{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Introduction\n\nIn the National Football League, pinning a team deep in their own territory or flipping the field with a deep, well-covered punt can lead to large momentum changes within a game.  The extra yardage gained by these punts can smother opposing offenses, help struggling defenses make important stops, as well as help your own offense get the ball back in better field position.  While these game changing punts are the amalgamation of several key components, such as protection and coverage schemes, an often overlooked component is the performance of the punter to set-up these plays.\n\nAs the \"quarterback\" of special team punts, the punter can read the return team's pre-snap coverage and determine where to punt the ball in order to minimize the opposing team's field position.  Ultimately, the value of interest is this resulting expected field position.  The location of the punt relative to everyone else on the field will determine the actions of the return team, i.e. fair catching the ball, attempting to return, or letting the ball bounce.  This decision by the return team then determines the final field position of the ball at the end of the punt play.\n\nImportantly, the analysis below is meant to evaluate the performance of the punter given only the pre-snap information available to the punter and the target landing location of the football in an attempt to decouple the play-specific performance of the other players.  Modelling the return team decision and the random bouncing of the football, each punt can be evaluated in terms of the resulting field position and compared to the expected value for an optimal punt location.","metadata":{}},{"cell_type":"markdown","source":"# 1. Data Processing\n\n### Where does the ball land?\n\nThe most consistent method to find the frame where the ball lands is to use the measured hang time and add this to the tracking data time from the punt frame.  For those plays where the 'punt' event frame is not confidently found, the punt frame can instead be found by adding the operation time from the time of the snap.  One method to check that the landing frame of the football has been consistently found is to the compare the position of the ball in this frame to the recorded kick length (and adding the original line of scrimmage).  When the ball is cleanly caught from the air, the two values are in clear agreement between, represented by the dashed line in the plot below.  However, the distance the ball bounces is not included when finding the landing frame of the ball, even though it is included in the kick length measurement.  As a result, balls bouncing backwards will have shorter kick lengths then the landing location, while conversely balls bouncing forward will appear longer.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom scipy.integrate import nquad\nimport seaborn as sns\n\nimport os\n\ndir_name = '/kaggle/input/nfl-big-data-bowl-2022/'\n\nplay_data = pd.read_csv(os.path.join(dir_name,'plays.csv'))\npff_data = pd.read_csv(os.path.join(dir_name,'PFFScoutingData.csv'))\n\ndiscard_results = ['Non-Special Teams Result']\npunt_plays = play_data[(play_data['specialTeamsPlayType']=='Punt')&(~play_data['specialTeamsResult'].isin(discard_results))]\npff_punt = pff_data.merge(punt_plays[['gameId','playId']])\n\n\n# 2018 Tracking Data\ntracking_game = pd.read_csv(os.path.join(dir_name,'tracking2018.csv'), chunksize=10**5)\nfootball_tracking = pd.concat((x.query(\"displayName == 'football'\") for x in tracking_game), ignore_index=True)\npunt_fbtrack_18 = football_tracking.merge(punt_plays[['gameId','playId']])\npunt_fbtrack_18['time'] = pd.to_datetime(punt_fbtrack_18['time'])\n\n# 2019 Tracking Data\ntracking_game = pd.read_csv(os.path.join(dir_name,'tracking2019.csv'), chunksize=10**5)\nfootball_tracking = pd.concat((x.query(\"displayName == 'football'\") for x in tracking_game), ignore_index=True)\npunt_fbtrack_19 = football_tracking.merge(punt_plays[['gameId','playId']])\npunt_fbtrack_19['time'] = pd.to_datetime(punt_fbtrack_19['time'])\n\n# 2020 Tracking Data\ntracking_game = pd.read_csv(os.path.join(dir_name,'tracking2020.csv'), chunksize=10**5)\nfootball_tracking = pd.concat((x.query(\"displayName == 'football'\") for x in tracking_game), ignore_index=True)\npunt_fbtrack_20 = football_tracking.merge(punt_plays[['gameId','playId']])\npunt_fbtrack_20['time'] = pd.to_datetime(punt_fbtrack_20['time'])\n\n# Combine all above tracking data\npunt_fbtrack = pd.concat([punt_fbtrack_18,punt_fbtrack_19,punt_fbtrack_20]) \\\n                .sort_values(by=['gameId','playId','frameId']) \\\n                .drop(columns=['o','dir','nflId','displayName','jerseyNumber','position','team']) \\\n                .reset_index(drop=True)\n\n# Add relevant PFF data\npunt_fbtrack = punt_fbtrack.merge(pff_punt[['gameId','playId','operationTime','hangTime','kickContactType']])\n\n\n\n# Select OOB and DEZ events and find all frames where the ball is out of bounds or past the endzones\n# Label these frames with a new event (if current event is None) to label that the ball has left the field of play.\noutbound_plays = punt_fbtrack[punt_fbtrack['kickContactType'].isin(['OOB','DEZ'])].copy()\noutbound_plays.loc[\n    ((outbound_plays['x']<10)|(outbound_plays['x']>110) # Touchback\n    |(outbound_plays['y']<0)|(outbound_plays['y']>160/3)) # Out of Bound\n    &(outbound_plays['event']=='None'), # Do not replace other events\n    'event'] = 'out_bound'\n\n# Find the frame that the ball goes out of bounds, or is caught/lands inbounds\noutbound_events = ['out_bound','out_of_bounds','punt_land','kick_received','punt_muffed','punt_received','fair_catch','touchback']\noutbound_plays_final = outbound_plays[outbound_plays['event'].isin(outbound_events)].drop_duplicates(subset=['gameId','playId'])\n\n\n\n# Add hangTime from PFF data to the frame where there is a punt or autoevent_kick event in the tracking data\n# This time gives a frame where the ball is recorded to land.\n# Most of the times, this frame will not have an event listed, but it will be very close to a large number of possible events\ninbound_plays = punt_fbtrack.merge(outbound_plays_final[['gameId','playId']], how='left', indicator=True).query('_merge == \"left_only\"').drop(columns=['_merge'])\n\n# Ball is not being tracked properly, so we cannot accurately find where the ball lands\nball_data_oob = inbound_plays[((inbound_plays['y']<0)|(inbound_plays['y']>160/3))&(inbound_plays['frameId']<30)].drop_duplicates(subset=['gameId','playId'])\ninbound_plays = inbound_plays.merge(ball_data_oob[['gameId','playId']],how='left',indicator = True).query('_merge == \"left_only\"').drop(columns=['_merge'])\n\n# Find punt frame and landing frame\npunt_frames = inbound_plays[inbound_plays['event'].isin(['punt','autoevent_kickoff'])].drop_duplicates(subset=['gameId','playId'],keep='last').copy()\npunt_frames['landFrame'] = punt_frames['frameId'] + np.ceil(10*punt_frames['hangTime'])\n\n# One play has the 'punt'-event frame clearly with significant error, so we exclude this event.\npunt_frames.drop(punt_frames[(punt_frames['gameId']==2018112507)&(punt_frames['playId']==560)].index,inplace=True)\n\n# Save the frames equal to the calculated landing frame\ninbound_plays = inbound_plays.merge(punt_frames[['gameId','playId','landFrame']], how='left')\nlf_from_ht = inbound_plays.query('frameId == landFrame')\n\n\n\n# Find the remaining in-bound plays\n# These plays are broken down into 3 sections\n# PFF has no hang time data, there was no 'Punt' event frame found, or the calculated landing frame is after the tracking data ends.\nib_remain = inbound_plays.merge(lf_from_ht[['gameId','playId']],how='left', indicator=True).query('_merge == \"left_only\"').drop(columns=['_merge']).copy()\nno_ht = ib_remain[ib_remain['hangTime'].isna()].copy()\nno_lt = ib_remain[(~ib_remain['hangTime'].isna())&(ib_remain['landFrame'].isna())].iloc[:,:-2].copy()\nno_frame = ib_remain[(~ib_remain['hangTime'].isna())&(~ib_remain['landFrame'].isna())].copy()\n\n# No \"punt\" frame event can be solved by adding operation time to the mix\nno_lt = no_lt.merge(pff_punt[['gameId','playId','operationTime']])\nsnap_frames = no_lt[no_lt['event']=='ball_snap'].copy()\nsnap_frames['landTime'] = snap_frames['frameId'] + np.floor(10*(no_lt['operationTime']+no_lt['hangTime']))\nno_lt = no_lt.merge(snap_frames[['gameId','playId','landTime']])\nno_lt_lf = no_lt.query('frameId == landTime')\n\n# When the landing frame is not present in the tracking data, we allow 3 or less frames of tolerance\n# This is likely similar error we have with bunching times in the tracking data\nno_frame_last = no_frame.drop_duplicates(subset=['gameId','playId'],keep='last')\nno_frame_last_keep = no_frame_last[no_frame_last['landFrame'] - no_frame_last['frameId']<=3]\nno_frame_last_remove = no_frame_last[no_frame_last['landFrame'] - no_frame_last['frameId']>3]\n\n# No hangtime data usually means the punt was blocked, which we will usually discard here.\n# However, some of these events were only deflected and still made it past the line of scrimmage\n# Last 2 bounce frames can be found manually watching the replay\n# gameId = 2018102111, playId = 3651, ball bounces on frame 38.\n# gameId = 2020100401, playId = 211,  ball bounces on frame 54.\nno_ht_def_punts = no_ht.merge(punt_plays[['gameId','playId','specialTeamsResult']]).query('specialTeamsResult != \"Blocked Punt\"')\nno_ht_lf = pd.concat([no_ht_def_punts[no_ht_def_punts['event'].isin(['fair_catch','punt_received','punt_land'])],\n    no_ht_def_punts[(no_ht_def_punts['gameId']==2018102111)&(no_ht_def_punts['playId']==3651)&(no_ht_def_punts['frameId']==38)],\n    no_ht_def_punts[(no_ht_def_punts['gameId']==2020100401)&(no_ht_def_punts['playId']==211)&(no_ht_def_punts['frameId']==54)]])\nno_ht_lf = no_ht_lf.merge(no_ht_def_punts[no_ht_def_punts['event']=='punt'][['gameId','playId','frameId']].rename({'frameId':'puntFrame'},axis=1))\nno_ht_lf['hangTime'] = (no_ht_lf['frameId'] - no_ht_lf['puntFrame'])/10  # Calculate time from punt frame to landing frame, this is definition of hangTime\n\n\n\n# Collect all landing frame data\nfb_landing = pd.concat([outbound_plays_final,lf_from_ht,no_lt_lf,no_frame_last_keep,no_ht_lf],ignore_index=True)[['gameId','playId','time','x','y','frameId','event','hangTime','playDirection']].rename({'x':'x_land','y':'y_land','frameId':'landFrame'},axis=1)\n\n\n# Compare the calculated landing location to the play data kick length + starting location\ncompare_xloc = punt_plays[['gameId','playId','absoluteYardlineNumber','kickLength']].merge(fb_landing[['gameId','playId','x_land','playDirection']]).merge(pff_punt[['gameId','playId','kickContactType']])\ncompare_xloc.loc[compare_xloc['playDirection']=='left', 'x_land'] = 120 - compare_xloc.loc[compare_xloc['playDirection']=='left', 'x_land']\ncompare_xloc.loc[compare_xloc['playDirection']=='left', 'absoluteYardlineNumber'] = 120 - compare_xloc.loc[compare_xloc['playDirection']=='left', 'absoluteYardlineNumber']\ncompare_xloc['x_tot'] = compare_xloc['absoluteYardlineNumber'] + compare_xloc['kickLength']\ncompare_xloc = compare_xloc[compare_xloc['kickContactType'].isin(['BF','CC', 'BB'])]\ncompare_xloc['dx'] = compare_xloc['x_tot'] - compare_xloc['x_land']\n\n\nlab_dict = {'BB':'Bounce Back', 'BF':'Bounce Forward', 'CC':'Clean Catch'}\n\nfig,ax = plt.subplots(figsize=(8,6))\nfor label,group_df in compare_xloc.groupby('kickContactType'):\n    ax.scatter(group_df['x_land'],group_df['x_tot'], label=lab_dict[label], alpha=0.5)\nax.legend(bbox_to_anchor=(.02, 1), loc='upper left')\nax.plot(np.linspace(0,120,2),np.linspace(0,120,2),'k--')\nplt.xlabel('Landing Location')\nplt.ylabel('Kick Length + Line of Scrimmage')\nplt.ylim(35,110)\nplt.xlim(20,120)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-06T16:39:44.459149Z","iopub.execute_input":"2022-01-06T16:39:44.459402Z","iopub.status.idle":"2022-01-06T16:41:27.448988Z","shell.execute_reply.started":"2022-01-06T16:39:44.459375Z","shell.execute_reply":"2022-01-06T16:41:27.448035Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Optimal 'League Average' Punting\n\nUltimately, the goal is to compare the performance of a punter to what a \"league-average\" punter would be capable.  While much of this averaging will be inherent in the modeling process below, one important aspect is to understand the distance a \"league-average\" punter is capable of kicking and how the hang time of each punt would change depending how far the punt travelled.\n\nThe complexity comes from not having access to information based on the vertical position of the football.  However, using kinematic equations, the distance the ball travels from the punter to the landing coordinates and the hang time are related to the velocity of the ball immediately following the punt by the following relation:\n\\begin{equation}\nv_\\text{punt}^2 = \\left(\\frac{\\Delta x}{T}\\right)^2 + \\left(g T\\right)^2,\n\\end{equation}\nwhere $\\Delta x$ is the distance the ball travels, $T$ is the hang time, and $g = 32.1$ ft$/s^2$ is the gravitational acceleration.  When considering a \"league average\" punter, the mean velocity of a kick can be found from the velocity distribution from all punts found using this formula.","metadata":{}},{"cell_type":"code","source":"###############################################################################################################\n# Find the snap frame for each play.  This allows some players to move to a final resting spot before the snap\n###############################################################################################################\n\nfb_use = punt_fbtrack.merge(fb_landing[['gameId','playId']]).copy()\nfb_use['dx'] = fb_use.groupby(['gameId','playId'])['x'].diff(-1)\n\n# Generally speaking, the 'ball_snap' event is very reliable, so if it is present in the data, always use that\nsnaps = fb_use[fb_use['event']=='ball_snap']\n\n# Otherwise, if the ball is already moving a significant distance (0.1 yards) in the first frame, then take the first frame as the snap frame\nno_snap = fb_use[~(fb_use['gameId'].astype(str) + fb_use['playId'].astype(str)) \\\n                 .isin((snaps['gameId'].astype(str) + snaps['playId'].astype(str)).values)]\nframe1_snap = no_snap[(no_snap['frameId']==1)&(np.abs(no_snap['dx'])>0.1)]\n\n#  For the remaining plays, take the first frame where the ball moves a significant amount as the snap frame\nremaining = fb_use.merge(no_snap[(no_snap['frameId']==1)&(np.abs(no_snap['dx'])<=0.1)][['gameId','playId']])\nearly_snap = remaining[np.abs(remaining['dx'])>0.1].groupby(['gameId','playId']).head(1)\n\nsnap_frames = pd.concat([snaps,frame1_snap,early_snap])[['gameId','playId','frameId']]\n\n\n\n############################################################\n# Pull the positional data for each player at the snap frame\n############################################################\n\nmax_frame = max(snap_frames['frameId'])\n\n# 2018 Tracking Data\ntracking_game = pd.read_csv(os.path.join(dir_name,'tracking2018.csv'), chunksize=10**5)\nfootball_tracking = pd.concat((x.query(\"frameId <= @max_frame\") for x in tracking_game), ignore_index=True)\npunt_fbtrack_18 = football_tracking.merge(snap_frames)\npunt_fbtrack_18['time'] = pd.to_datetime(punt_fbtrack_18['time'])\n\n# 2019 Tracking Data\ntracking_game = pd.read_csv(os.path.join(dir_name,'tracking2019.csv'), chunksize=10**5)\nfootball_tracking = pd.concat((x.query(\"frameId <= @max_frame\") for x in tracking_game), ignore_index=True)\npunt_fbtrack_19 = football_tracking.merge(snap_frames)\npunt_fbtrack_19['time'] = pd.to_datetime(punt_fbtrack_19['time'])\n\n# 2020 Tracking Data\ntracking_game = pd.read_csv(os.path.join(dir_name,'tracking2020.csv'), chunksize=10**5)\nfootball_tracking = pd.concat((x.query(\"frameId <= @max_frame\") for x in tracking_game), ignore_index=True)\npunt_fbtrack_20 = football_tracking.merge(snap_frames)\npunt_fbtrack_20['time'] = pd.to_datetime(punt_fbtrack_20['time'])\n\n# Combine all above tracking data\ntrack_snap = pd.concat([punt_fbtrack_18,punt_fbtrack_19,punt_fbtrack_20]) \\\n                .sort_values(by=['gameId','playId','frameId']) \\\n                .drop(columns=['s','a','dis','event','frameId','displayName','jerseyNumber']) \\\n                .reset_index(drop=True)\n\n\n# Flip tracking data so that all punts are 'right'-directional, i.e. the punt moves towards larger x-positions.\n\ndef flip_play_direction(df):\n    df_flipped = df.copy()\n    df_flipped.loc[df_flipped['playDirection']=='left', 'x'] = 120 - df_flipped.loc[df_flipped['playDirection']=='left', 'x']\n    df_flipped.loc[df_flipped['playDirection']=='left', 'y'] = (160/3 - df_flipped.loc[df_flipped['playDirection']=='left', 'y']).round(2)\n    \n    # When orientation of players is included, flipping field is equal to a 180 degree rotation.\n    # Mod by 360 degrees to keep values in the range [0,360)\n    orientation_cols = [col for col in df_flipped.columns if col in ['o','dir']]\n    for col in orientation_cols:\n        df_flipped[col] = (180*(df_flipped['playDirection']=='left') + df_flipped[col])%360\n    \n    return df_flipped\n\nfbland_f = flip_play_direction(fb_landing.rename({'x_land':'x', 'y_land':'y'},axis=1))\nsnap_track_f = flip_play_direction(track_snap)\n\nfb_snap = snap_track_f[snap_track_f['team']=='football'].copy() # Football\nplayer_snap = snap_track_f[snap_track_f['team'].isin(['home','away'])] # Players at snap frame\n\n# Every play has a designated punters (either 'P' or 'K')\npunter = player_snap[player_snap['position'].isin(['P','K'])].sort_values('x',ascending=False).groupby(['gameId','playId']).last().reset_index()\n\n\n# Calculate punter's kick velocity\nkick_length = punter[['gameId','playId','x','y','nflId']].rename({'x':'x_punt', 'y':'y_punt'},axis=1).merge(fbland_f[['gameId','playId','x','y','hangTime']].rename({'x':'x_land', 'y':'y_land'},axis=1))\nkick_length['dx'] = kick_length['x_land'] - kick_length['x_punt']\nkick_length['dy'] = kick_length['y_land'] - kick_length['y_punt']\nkick_length['dist'] = np.sqrt(kick_length['dx']**2 + kick_length['dy']**2)\n\ndef calc_const(dist, ht):\n    g = 32.1/3\n    v_const = (dist/ht)**2 + (g*ht/2)**2\n    return np.sqrt(v_const)\n\nkick_length['v_punt'] = calc_const(kick_length['dist'],kick_length['hangTime'])\nvc_avg = kick_length['v_punt'].describe()['mean']\n\n\ndef calc_time(del_x, v_const):\n    g = 32.1/3\n    p1 = 2/g**2*v_const**2\n    p2 = 2/g**2*np.sqrt(v_const**4 - g**2 * del_x**2)\n    return [np.sqrt(p1-p2),np.sqrt(p1+p2)]\n\n\nfig, ax = plt.subplots(1,2, figsize=(12,6))\nsns.histplot(data=kick_length, x='v_punt', stat='density', ax=ax[0])\nax[0].axvline(x=vc_avg, label=f'Mean = {vc_avg:.1f} yd/s',c='k', ls='--')\nax[0].legend()\nax[0].set_xlim(15,32.5)\n\ndist_array = np.arange(0,3*vc_avg**2/32.1,0.1)\nht_m, ht_p = calc_time(dist_array, vc_avg)\nax[1].scatter(dist_array, ht_p, alpha=0.5)\nax[1].scatter(dist_array, ht_m, alpha=0.5)\nax[1].plot(np.array([0,70]), np.sqrt(2/(32.1/3)**2*vc_avg**2)*np.array([1,1]), 'k--')\nax[1].text(24, 4.45, 'High launch angle', fontsize=12)\nax[1].text(24, 2, 'Low launch angle', fontsize=12)\nax[1].text(42,0.1, 'velocity = 27.2 yd/s', fontsize=12)\nax[1].set_xlabel('Distance (yds)')\nax[1].set_ylabel('Hang Time (s)')\nax[1].set_xlim(0,70)\nax[1].set_ylim(0,5.2)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-06T16:41:27.451328Z","iopub.execute_input":"2022-01-06T16:41:27.451883Z","iopub.status.idle":"2022-01-06T16:42:39.765802Z","shell.execute_reply.started":"2022-01-06T16:41:27.451838Z","shell.execute_reply":"2022-01-06T16:42:39.764936Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using this kinematic relation, given a punt velocity and the punter's original location, any position on the field is a potential landing location given the condition\n\\begin{equation}\nT = \\frac{2}{g^2} \\left[v_\\text{punt}^2 + (v_\\text{punt}^4 - g^2 \\Delta x^2)^{1/2} \\right] \\;\\; \\rightarrow \\;\\; \\Delta x \\leq \\frac{v_\\text{punt}^2}{g}.\n\\end{equation}\nThis hang time is always taken as the higher launch angle solution, even if the real punt is calculated at the shorter solution.","metadata":{}},{"cell_type":"markdown","source":"# 2. Fair Catch, Return, or Not Fielded?\n\nIn order to model the return team's decision making process on whether or not to field the football, a 2D Convolution Neural Network is used. The general model structure is very similar to the work used by The Zoo team in the 2020 Big Data Bowl (write up here: https://www.kaggle.com/c/nfl-big-data-bowl-2020/discussion/119400) modeling the length of a rushing play.  The two dimensions of the CNN represent punt team vs return team data, several features encoding the relative position data of the team.  In their work, the player/location of central interest was the running back, which acts as an attractor for all defensive players.  However for a punt, some players on the return team are attempting to pressure the punter, while others are moving towards the landing position of the ball.","metadata":{}},{"cell_type":"code","source":"inbound_land = fb_landing[(fb_landing['x_land']<110)&(fb_landing['x_land']>10)&(fb_landing['y_land']>0)&(fb_landing['y_land']<160/3)].copy()\ninbound_land = inbound_land.merge(punt_plays[['gameId','playId','specialTeamsResult']]).merge(pff_punt[['gameId','playId','kickContactType']]).sort_values(['gameId','playId'])\n\n\n############################################################################\n# Assign each punt to be either a bounce or a catch (For binary classifier)\n############################################################################\n\n# For the plays with the ball landing in-bounds, we need to categorize whether the receiving team catches (or attempts to catch) the punt or not.\n# This can be done looking at the Contact Type recorded for each punt.  Some require more nuance.\n\n# All punts of a particular contact type are bounces or catches\ncontact_to_land = {\n    'BB':'Bounce',\n    'BC':'Catch',\n    'BF':'Bounce',\n    'BOG':'Bounce',\n    'CC':'Catch',\n    'CFFG':'Bounce',\n    'DEZ':'Bounce',\n    'KTB':'Bounce',\n    'KTC':'Bounce',\n    'KTF':'Bounce',\n}\n\n# The remaining punts (ICC, MBC, MBDR, OOB, NaN) need to be looked at individually.\n# OOB with the event as out_of_bounds should be listed as out_plays, so we must remove.\n# For OOB and NaN events, the events are all filled and correspond to the following catch/bounce result:\nevent_to_land = {\n    'kick_received':'Catch',\n    'punt_received':'Catch',\n    'punt_land':'Bounce',\n    'fair_catch':'Catch',\n    'fumble_defense_recovered':'Catch',\n    'punt_muffed':'Catch',\n    'fumble':'Catch' # For MBDR event\n}\n\n# Look at MBC punts:\n# 2018112504 - 2312 -> 'Catch' (Found from the replay highlights)\n# All others are 'Bounce'\n\n# Look at ICC punts:\n# specialTeamsResult = Fair Catch -> 'Catch'\n\n\n# Remaining events, make event table look at [2] in event table\nevent_to_land_remain = {\n    'punt_received':'Catch',\n    'first_contact':'Catch',\n    'punt_muffed':'Catch',\n    'punt_received':'Catch',\n    'punt_downed':'Bounce',\n    'punt_land':'Bounce',\n    'fair_catch':'Catch',\n    'kick_received':'Catch',\n    'fumble':'Catch'\n}\n\n\n# Create ordered list of events (exluding those in a particular list) for each play\ndef event_table(df, exc_list):\n    df_trim = df[(~df['event'].isin(exc_list))&(df['event']!='None')].copy()\n    df_trim['event_num'] = df_trim.groupby(['gameId','playId']).cumcount().values\n    return pd.pivot_table(df_trim, values='event', index=['gameId','playId'], columns='event_num', aggfunc=lambda x: ''.join(x))\n\n\n# Remove plays that are called out of bounds, even if the landing frame claims they are in bound\n# This is mostly for consistency.  Either this landing frame is found incorrectly, or the refs whistled the play dead, so the return team may play differently\ninbound_land_remoob = inbound_land.drop(inbound_land[((inbound_land['kickContactType']=='OOB')|(inbound_land['kickContactType'].isna()))&(inbound_land['event']=='out_of_bounds')].index)\n\n# Using the above dictionaries, assign Catch/Bounce labels to each punt (without needing event table)\ninbound_land_remoob['Type'] = inbound_land_remoob['kickContactType'].map(contact_to_land)\ninbound_land_remoob.loc[inbound_land_remoob['Type'].isna(),'Type'] = inbound_land_remoob[(inbound_land_remoob['kickContactType'].isin(['OOB','MBDR']))|(inbound_land_remoob['kickContactType'].isna())]['event'].map(event_to_land)\ninbound_land_remoob.loc[inbound_land_remoob['kickContactType']=='MBC','Type'] = 'Bounce'\ninbound_land_remoob.loc[(inbound_land_remoob['gameId']==2018112504)&(inbound_land_remoob['playId']==2312),'Type'] = 'Catch'\ninbound_land_remoob.loc[(inbound_land_remoob['kickContactType']=='ICC')&(inbound_land_remoob['specialTeamsResult']=='Fair Catch'),'Type'] = 'Catch'\n\n# Use event table and the last dictionary to assign Catch/Bounce, then merge the results to have the full list of punt results.\nremaining_ibland = inbound_land_remoob[inbound_land_remoob['Type'].isna()]\nfull_event_punts_table = event_table(punt_fbtrack.merge(remaining_ibland[['gameId','playId']]),['line_set']).reset_index([0,1]).merge(remaining_ibland[['gameId','playId','specialTeamsResult','kickContactType']])\nfull_event_punts_table['Type_update'] = full_event_punts_table[2].map(event_to_land_remain)\n\ninbound_land_remoob = inbound_land_remoob.merge(full_event_punts_table[['gameId','playId','Type_update']], how='left')\ninbound_land_remoob.loc[inbound_land_remoob['Type'].isna(), 'Type'] = inbound_land_remoob.loc[inbound_land_remoob['Type'].isna(), 'Type_update']\ninbound_land_remoob.drop(columns=['Type_update', 'specialTeamsResult','kickContactType','event'],inplace=True)\n\nbounce_catch = inbound_land_remoob.sort_values(['gameId','playId'])","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-06T16:42:39.767252Z","iopub.execute_input":"2022-01-06T16:42:39.767551Z","iopub.status.idle":"2022-01-06T16:42:39.903902Z","shell.execute_reply.started":"2022-01-06T16:42:39.767511Z","shell.execute_reply":"2022-01-06T16:42:39.90307Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Separate players based on whether they are on the same team as the punter or not\n# Note, we remove the punter from the punt team right away.\ntemp = player_snap.merge(punter[['gameId','playId','nflId','team']],on=['gameId','playId'],suffixes=('','_punt'))\npunt_team = temp[(temp['team']==temp['team_punt'])&(temp['nflId']!=temp['nflId_punt'])].copy()\nrec_team = temp[temp['team']!=temp['team_punt']].copy()\n\n# For each play, get a list of x and y positions for each player at the snap for both punt and return team\npunt_team['num'] = punt_team.groupby(['gameId','playId']).cumcount()\npunt_team_pos = pd.pivot_table(punt_team,values=['x','y'],index=['gameId','playId'],columns='num')\npunt_team_pos.columns  = [f'x_p{i}' for i in range(10)] + [f'y_p{i}' for i in range(10)]\npunt_team_pos = punt_team_pos[[punt_team_pos.columns[i] for i in np.array([[j, j+10] for j in range(10)]).reshape((20,))]].reset_index()\n\nrec_team['num'] = rec_team.groupby(['gameId','playId']).cumcount()\nrec_team_pos = pd.pivot_table(rec_team,values=['x','y'],index=['gameId','playId'],columns='num')\nrec_team_pos.columns  = [f'x_r{i}' for i in range(11)] + [f'y_r{i}' for i in range(11)]\nrec_team_pos = rec_team_pos[[rec_team_pos.columns[i] for i in np.array([[j, j+11] for j in range(11)]).reshape((22,))]].reset_index()\n\n# Get location of football at the snap.  Note, we do not use the tracking data for the x-value of the football\n# This is because sometimes the snap frame is slightly delayed, and we really want to know what the original line of scrimmage is\nfootball_loc = fb_snap[['gameId','playId','x','y','playDirection']].merge(punt_plays[['gameId','playId','absoluteYardlineNumber']]) \\\n                    .rename({'absoluteYardlineNumber':'x_yl'},axis=1)\nfootball_loc['x_yl'] += (football_loc['playDirection']=='left')*(120-2*football_loc['x_yl'])\nfb_loc = football_loc[['gameId','playId','x_yl','y']].rename({'x_yl':'x_fb','y':'y_fb'},axis=1)\n\n\n\n# Combine all the individual play data needed to train the fielding decision model\nall_pos = punter[['gameId','playId','nflId','x','y']].rename({'nflId':'nflId_punt','x':'x_punt','y':'y_punt'},axis=1)\\\n    .merge(punt_team_pos).merge(rec_team_pos)\\\n    .merge(fb_loc).merge(fbland_f[['gameId','playId','x','y','hangTime']].rename({'x':'x_land','y':'y_land'},axis=1)) \\\n    .merge(bounce_catch[['gameId','playId','Type']]).merge(punt_plays[['gameId','playId','specialTeamsResult']])\n\n\n\n# This final dataframe contains all the information used to train the CNNs.  The individual training is saved for separate notebooks linked above.","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-06T16:42:39.905577Z","iopub.execute_input":"2022-01-06T16:42:39.905885Z","iopub.status.idle":"2022-01-06T16:42:40.131075Z","shell.execute_reply.started":"2022-01-06T16:42:39.905854Z","shell.execute_reply":"2022-01-06T16:42:40.130398Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The nine features used in the model are simply the relative x and y positions of the players or locations of interest, in addition to the hang time of the punt.  Note, the punter is excluded from the punt team, giving the result 11 x 10 x 9 tensor with which to train.\n\nFeatures:\n* (Punt Team - Return Team) X, Y\n* (Punter - Return Team) X, Y\n* (Landing Location - Return Team) X, Y\n* Presnap Football Location X, Y\n* Hang time\n\nWith these features, two separate networks are trained: One trained for the binary classification of whether the punt will be fielded or not; the other trained on classifying the three outcomes of a punt being returned, fair caught, or not fielded, which are then combined taking an ensemble average.  The goal of training the models in this way is to encourage the distinction between fielding and not fielding the punt.  Explicitly, the binary model outputs the probability of the punt not being fielded, while the three-class model outputs the probabilities of being a fair catch, not fielded, and a return.\n\nthe output of these two models are combined to find the following probabilities:\n\\begin{align}\n\\text{Not Fielded}&: P(\\text{NF}) = \\frac{1}{2}\\left(P_{\\text{binary}} + P_{\\text{3-out}}(\\text{NF})\\right)\\\\\n\\text{Fair Catch}&: P(\\text{FC}) = \\left(1 - P(\\text{NF})\\right) \\times \\left(\\frac{P_{\\text{3-out}}(\\text{FC})}{P_{\\text{3-out}}(\\text{FC}) + P_{\\text{3-out}}(\\text{R})}\\right) \\\\\n\\text{Return}&: P(\\text{R}) = \\left(1 - P(\\text{NF})\\right) \\times \\left(\\frac{P_{\\text{3-out}}(\\text{R})}{P_{\\text{3-out}}(\\text{FC}) + P_{\\text{3-out}}(\\text{R})}\\right)\n\\end{align}\n\nThe above ensemble can be used to predict the probabilities for each outcome at each location on the field, making sure to modify the hang time for each proposed punt.  In the example plot below, the punt team is backed up deep in their own territory, with both gunners double covered.  The return team is attempting to take advantage of the already good field position and gain even more yardage with a return.","metadata":{"_kg_hide-output":true}},{"cell_type":"code","source":"g = 2020121313 # 2018100703\np = 2577 # 4645\nfbloc_ensprob = pd.read_csv(f'/kaggle/input/exampleensemble/ExampleEnsembleProb_{g}_{p}.csv')\nens_pos_ex1 = all_pos[(all_pos['gameId']==g)&(all_pos['playId']==p)]\n\nfb_x = ens_pos_ex1.iloc[:,[47,49]]\nfb_y = ens_pos_ex1.iloc[:,[48,50]]\n\npunter_loc = ens_pos_ex1.iloc[:,3:5]\n\npunt_team_x = ens_pos_ex1.iloc[:,5:24:2]\npunt_team_y = ens_pos_ex1.iloc[:,6:25:2]\n\nret_team_x = ens_pos_ex1.iloc[:,25:46:2]\nret_team_y = ens_pos_ex1.iloc[:,26:47:2]\n\n\n# Hash mark variables\nhash_up_mid = 160/3 - 70.75/3\nhash_down_mid = 70.75/3\nhash_length = 2/3\nx_hash_lines = np.arange(11,110,1)\ny_hash_lines_up = [hash_up_mid-hash_length/2, hash_up_mid+hash_length/2]\ny_hash_lines_down = [hash_down_mid-hash_length/2, hash_down_mid+hash_length/2]\n\ngs_kw = dict(width_ratios=[8, 1], height_ratios=[8, 1])\nfig, ax = plt.subplot_mosaic([['left', 'ul'],['left','ll']], gridspec_kw=gs_kw, figsize=(15, 6))\n\n\n# Plot model output\nax['left'].scatter(fbloc_ensprob['x_land'],fbloc_ensprob['y_land'],c=fbloc_ensprob[['R','NF','FC']].values)\n\n# Draw hash marks and 5-yard intervals\nhash_alpha = 0.2\nax['left'].vlines(x=x_hash_lines,ymin=y_hash_lines_up[0],ymax=y_hash_lines_up[1],colors='k',alpha=hash_alpha) # Upper hash\nax['left'].vlines(x=x_hash_lines,ymin=y_hash_lines_down[0],ymax=y_hash_lines_down[1],colors='k',alpha=hash_alpha) # Lower hash\nax['left'].vlines(x=np.arange(15,110,5),ymin=0,ymax=54,colors='k',alpha=hash_alpha) # 5-yard intervals\nax['left'].vlines(x=[12,108],ymin=80/3-0.5,ymax=80/3+0.5,colors='k',alpha=hash_alpha) # 2-point conversion\nax['left'].vlines(x=[10,110],ymin=0,ymax=55,colors='k') # Endzones\n#plt.vlines(x=yardline,ymin=0,ymax=55,colors='k',linestyles='dashed') # Line of Scrimmage\n\n# Add text yard marker\nax['left'].text(20,12,'10',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(30,12,'20',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(40,12,'30',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(50,12,'40',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(60,12,'50',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(70,12,'40',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(80,12,'30',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(90,12,'20',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\nax['left'].text(100,12,'10',horizontalalignment='center',verticalalignment='top',fontsize=30,alpha=hash_alpha)\n\nax['left'].text(20-1,160/3-12,'10',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(30-1,160/3-12,'20',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(40-0.5,160/3-12,'30',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(50-0.5,160/3-12,'40',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(60-0.8,160/3-12,'50',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(70-0.5,160/3-12,'40',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(80-0.5,160/3-12,'30',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(90-1,160/3-12,'20',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\nax['left'].text(100-1,160/3-12,'10',horizontalalignment='center',verticalalignment='bottom',fontsize=30,alpha=hash_alpha,rotation='180')\n\n# Plot each team's players\nax['left'].scatter(pd.concat([punt_team_x,punter_loc['x_punt']],axis=1).iloc[0],pd.concat([punt_team_y,punter_loc['y_punt']],axis=1).iloc[0],c='k',marker='x',s=50)\nax['left'].scatter(ret_team_x.iloc[0],ret_team_y.iloc[0],s=50,facecolors='none',edgecolors='k')\n\nax['left'].set_xlim(0,120)\nax['left'].set_ylim(0,160/3)\nax['left'].set_aspect(1)\n\n\n\ntest = np.arange(0,1,0.002)\ntest = pd.DataFrame([[x,y] for x in test for y in test],columns=['r','g'])\ntest['b'] = 1 - test['r'] - test['g']\ntest = test[test['b']>=0].round(2)\n\nbasis = np.array([[0.0, 1.0], [-1.5/np.sqrt(3), -0.5],[1.5/np.sqrt(3), -0.5]])\nxy = np.dot(test.values,basis)\ntest['x'] = xy[:,0]\ntest['y'] = xy[:,1]\n\noffset = 0.5\nfontsize = 12\n\nax['ul'].scatter(test['x'],test['y'], c=test[['r','g','b']].values)\nax['ul'].text(basis[0,0]*(1+offset), basis[0,1]*(0.8+offset), 'P(R)', horizontalalignment='center',\n            verticalalignment='center', fontsize=fontsize)\nax['ul'].text(basis[1,0]*(1+offset), basis[1,1]*(1+offset), 'P(NF)', horizontalalignment='center',\n            verticalalignment='center', fontsize=fontsize)\nax['ul'].text(basis[2,0]*(1+offset), basis[2,1]*(1+offset), 'P(FC)', horizontalalignment='center',\n            verticalalignment='center', fontsize=fontsize)\n\nax['ll'].text(0.4, 0.5, f'Game ID: {g}', transform=ax['ll'].transAxes,ha=\"center\", va=\"center\", fontsize=14, color=\"k\")\nax['ll'].text(0.4, 0.0, f'Play ID: {p}', transform=ax['ll'].transAxes,ha=\"center\", va=\"center\", fontsize=14, color=\"k\")\n\nax['ul'].set_frame_on(False)\nax['ul'].set_xticks(())\nax['ul'].set_yticks(())\nax['ul'].set_aspect(1)\nax['ll'].set_frame_on(False)\nax['ll'].set_xticks(())\nax['ll'].set_yticks(())\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-06T16:42:40.132373Z","iopub.execute_input":"2022-01-06T16:42:40.132716Z","iopub.status.idle":"2022-01-06T16:42:43.18358Z","shell.execute_reply.started":"2022-01-06T16:42:40.132689Z","shell.execute_reply":"2022-01-06T16:42:43.182976Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Return Length\n\nThe model for return length is similar in structure to the previous model, using a 2d CNN over punt vs return team player location data.  One significant difference however is removing the presnap location of the ball from the features.  Intuitively, assuming a punt is caught and returned, only relative positions should have a significant effect on the return length, outside of touchdown returns or safeties.\n\nFurthermore, there is no additional feature indicating which player is the returner.  The main reason for this is that who will return the ball cannot always be known a priori.  There are a significant number of events with multiple members of the return team back deep in the field, and which member will return the ball may be dependent on where the ball is placed.  Similarly, there are laterals and longer throws that occasionally arise from the return team.  All of these add complexity to the return length (given only presnap information), many of which may be revealed in the formation.  Separating out a specific returner may miss out on these features.\n\nAs mentioned before, there is bound to be significant loss in any particular model using only the limited information available to the punter pre-snap.  Information such as blown lane assignments by the punt team can develop early within the punt and indicate potential large returns.  At the same time, these general mistakes are not (entirely) the fault of the punter and cannot be attributed as such.  Instead, the punt should evaluated assuming minimal mistakes, such as no missed tackles or muffed punts.  While there is certainly merit in some punts being more difficult to catch than others, leading to more mistakes by the returner, such considerations are beyond the scope of the present work.","metadata":{}},{"cell_type":"code","source":"outcome_land_data = fb_landing[['gameId','playId','landFrame']].merge(bounce_catch[['gameId','playId','Type']]).merge(play_data[['gameId','playId','specialTeamsResult']]).copy()\n\n# Find punts that were returned, even if they bounced first\nbounce_nf = ['Downed','Out of Bounds','Touchback']\npunt_return = ['Return','Out of Bounds','Blocked Punt','Muffed']\noutcome_land_data['Label'] = outcome_land_data['Type']\noutcome_land_data.loc[(outcome_land_data['Type']=='Bounce')&(~outcome_land_data['specialTeamsResult'].isin(bounce_nf)),'Label'] = 'Catch'\noutcome_land_data.loc[(outcome_land_data['Label']=='Catch')&(outcome_land_data['specialTeamsResult'].isin(punt_return)),'Label'] = 'Return'\nreturnable = outcome_land_data[(outcome_land_data['Label']=='Return')&(~outcome_land_data['specialTeamsResult'].isin(['Muffed','Downed','Fair Catch']))].drop(columns=['Label'])\n\n# Look at tracking data from frames after the ball has landed\nafter_land = punt_fbtrack.merge(returnable).query('frameId >= landFrame')\nframes_wevents = after_land.query('event != \"None\"')\nexcl_events=['punt_received','punt_land','kick_received','fair_catch','penalty_flag','lateral','pass_outcome_caught']\nframes_wevents = frames_wevents[~frames_wevents['event'].isin(excl_events)]\nframes_wevents['count'] = frames_wevents.groupby(['gameId','playId']).cumcount()\n\n# After exluding a large number of events, nearly always the first event after the landing frame is first contact.\n# The two fumble recovery events are:\n# 2018 game > go to first contact\n# 2019 game > go to recover\n# Note there is also an infamous touchback/safety/return punt caused by unsportsmanlike conduct penalty, which we remove\nreturn_values = frames_wevents.loc[(frames_wevents['count']==0)&(frames_wevents['event']!='touchback'),['gameId','playId','frameId','x','y','event']].copy()\nreturn_values.loc[return_values['event']=='fumble_defense_recovered',['frameId','x','y']] = \\\n       frames_wevents.merge(frames_wevents[(frames_wevents['event']=='fumble_defense_recovered')\n                                         &(frames_wevents['count']==0)][['gameId','playId']]).query('event==\"first_contact\"')[['frameId','x','y']].values\n\n# Given these first contact frames, calculate the yards gained from the landing location to first contact\nreturn_data = return_values.merge(after_land[['gameId','playId','playDirection']].drop_duplicates()).reset_index(drop=True).sort_values(['gameId','playId'])\nreturn_dist = flip_play_direction(return_data).merge(fbland_f[['gameId','playId','x','y']], on=['gameId','playId'],suffixes=['','_land']).merge(punt_plays[['gameId','playId','kickReturnYardage']]).copy()\nreturn_dist['YardsToContact'] = return_dist['x_land'] - return_dist['x']\n\nsns.displot(data=return_dist, x='YardsToContact', y='kickReturnYardage')\nplt.xlim(-20,100)\nplt.ylim(-20,100)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-06T17:15:53.800816Z","iopub.execute_input":"2022-01-06T17:15:53.801138Z","iopub.status.idle":"2022-01-06T17:15:54.336298Z","shell.execute_reply.started":"2022-01-06T17:15:53.801103Z","shell.execute_reply":"2022-01-06T17:15:54.335397Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To this end, assuming minimal mistakes by the punt team, the return length will instead be replaced by the yards gained from the landing location to first contact.  Only in rare occasions is there a significant deviation between the yards to first contact and the total return yards.  Large deviations between these two metrics are generally due to either a penalty, missed tackle, or blown coverage.\n","metadata":{}},{"cell_type":"markdown","source":"# 4. Bouncing Statistics\n\nIn much the same way a coin flip is considered random, the shape of the football enables the bouncing of the football to be largely considered random.  Nevertheless, there is an element of skill in the random bouncing of the ball.  In particular, punting the ball so that the gunners are able to stop unfavorable bounces can enable many punts near the endzone to be downed inside the five yardline.  Furthermore, the spin assocaited with the ball imparted by the handedness of the kicker can have a significant effect on biasing the bouncing direction and magnitude.\n\nThe model assumes that the x and y distances of the bounce are drawn from a bivariate normal distribution, but both means and the elements of the covariance matrix are functions of some features of the punt,\n\\begin{equation}\nX \\sim \\mathcal{N}(\\mu(\\text{features}), \\Sigma(\\text{features})).\n\\end{equation}\nThe direction the punt is travelling in the air should already have some forward momentum and should preferentially favor bouncing in that direction.  (Similarly, the spin orientation of the ball will bias the perpendicular direction, but this data is not available.)  It seems natural then to attempt to reduce the covariance of the two directions as much as possible by working in the coordinates given by this punt direction and its perpendicular.","metadata":{}},{"cell_type":"code","source":"bounce_punts = bounce_catch.query('Type == \"Bounce\"')\nframes_bounce = punt_fbtrack.merge(fb_landing[['gameId', 'playId','landFrame']]).query('frameId >= landFrame').merge(bounce_punts[['gameId','playId']])\n\nbounce_play_data = punt_plays[['gameId','playId','kickLength','penaltyYards','specialTeamsResult']].merge(bounce_punts[['gameId','playId']])\n\n# Frames where there is an easy to find event indicating some interaction with the ball after the bounce\nground_event_list = ['punt_downed','out_of_bounds','touchback','punt_received','punt_muffed','fair_catch','kick_received']\nfinal_event = frames_bounce[frames_bounce['event'].isin(ground_event_list)].drop_duplicates(subset=['gameId','playId'])\nremain = frames_bounce.merge(final_event[['gameId','playId']],indicator=True,how='outer').query('_merge == \"left_only\"').drop(columns=['_merge'])\n\nresult_to_event = {'Downed':'punt_downed','Fair Catch':'fair_catch','Muffed':'punt_muffed','Out of Bounds':'out_of_bounds','Touchback':'touchback','Return':'Return'}\nfe_merged = final_event.merge(bounce_play_data)\nfe_merged['event_check'] = fe_merged['specialTeamsResult'].map(result_to_event)\n\n\n# (5) Result = Muffed > Event = out_of_bounds ------ See below\n# gameId = 2018092000,playId = 801 : event = fumble_defense_recovered\n# gameId = 2018093005,playId = 3731 : frame = 88\n# gameId = 2018102104,playId = 789 : ?\n# gameId = 2018122309,playId = 3195 : ?\n# gameId = 2020101112,playId = 1389 : ?\n# Drop these events entirely.  Generally the ball does not bounce for long before being muffed\n\n# (1) Result = Downed > Event = fair_catch --------- 'Fake' fair catch -> skip to downed event\n# (5) Result = Downed > Event = punt_received ------ 2019: 'Fake' received (3) -> skip to downed event. 2020: Correct values (2)\n\n# All of the rest of the mismatched events are correct.\ncorrect = fe_merged.query('(event_check == event)|(event_check == \"Return\")')\ncheck_events = fe_merged.query('(event_check != event)&(event_check != \"Return\")')\n\ndropped_events = check_events.query('(specialTeamsResult == \"Muffed\")&(event == \"out_of_bounds\")')\nreeval_events = check_events[(check_events['specialTeamsResult']=='Downed')&((check_events['event']==\"fair_catch\")|((check_events['event']==\"punt_received\")&(check_events['gameId']//10**6 == 2019)))]\n\nchecked = check_events.merge(dropped_events,how='outer',indicator=True).query('_merge == \"left_only\"').drop(columns=['_merge']) \\\n    .merge(reeval_events,how='outer',indicator=True).query('_merge == \"left_only\"').drop(columns=['_merge'])\ncorrect = pd.concat([correct,checked,after_land.merge(reeval_events[['gameId','playId']]).query('event == \"punt_downed\"')])\n\nend_bounce = correct[['gameId','playId','x','y','frameId']].copy()\nfb_bounce_data = end_bounce.merge(fb_landing,on=['gameId','playId'])\nfb_bounce_data = flip_play_direction(fb_bounce_data).rename({'x':'x_final', 'y':'y_final'},axis=1)[['gameId','playId','x_land','y_land','x_final','y_final','landFrame','frameId','hangTime','playDirection']]\n\n# Gather data to find the bounce length in each direction\nfb_locations = fb_loc.merge(punter[['gameId','playId','x','y']]).merge(fb_bounce_data).rename({'x':'x_punt', 'y':'y_punt'}, axis=1).copy()\nfb_locations['dx_air'] = fb_locations['x_land']-fb_locations['x_punt']\nfb_locations['dy_air'] = fb_locations['y_land']-fb_locations['y_punt']\nfb_locations['dist_air'] = np.sqrt(fb_locations['dx_air']**2 + fb_locations['dy_air']**2).round(2)\nfb_locations['v_c'] = calc_const(fb_locations['dist_air'],fb_locations['hangTime'])\nfb_locations['v_field_ratio'] = fb_locations['dist_air']/fb_locations['hangTime'] / fb_locations['v_c']\n\nfb_locations['dx_bounce'] = fb_locations['x_final']-fb_locations['x_land']\nfb_locations['dy_bounce'] = fb_locations['y_final']-fb_locations['y_land']\nfb_locations['dist_bounce'] = np.sqrt(fb_locations['dx_bounce']**2 + fb_locations['dy_bounce']**2).round(2)\n\n# Coordinate rotation from x-y into the direction of the punt and its perpendicular axis\nfb_locations['cos_air_bounce'] = (fb_locations['dx_air']*fb_locations['dx_bounce'] + fb_locations['dy_air']*fb_locations['dy_bounce'])/(fb_locations['dist_air'])#*fb_locations['dist_bounce'])\nfb_locations['sin_air_bounce'] = (-fb_locations['dy_air']*fb_locations['dx_bounce'] + fb_locations['dx_air']*fb_locations['dy_bounce'])/(fb_locations['dist_air'])#*fb_locations['dist_bounce'])\n\nfb_locations = fb_locations.fillna(0).drop(columns=['x_land','y_land','x_final','y_final','landFrame','frameId','playDirection'])","metadata":{"execution":{"iopub.status.busy":"2022-01-06T09:08:45.47901Z","iopub.execute_input":"2022-01-06T09:08:45.479326Z","iopub.status.idle":"2022-01-06T09:08:45.946545Z","shell.execute_reply.started":"2022-01-06T09:08:45.479293Z","shell.execute_reply":"2022-01-06T09:08:45.945851Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This model is implemented using a Mixture Density Network.  The network will output the mean and standard deviation in each direction and correlation between them which minimizes the negative log likelihood.  The features for the model are given as:\n* Original line of scrimmage\n* Y-distance travelled in the air (punter_y - land_y)\n* Hang time\n* Magnitude of the (3d) velocity of the punt, as found in Section 1 Eq 1\n* Ratio of in-field velocity to total velocity\n\nPhysically, the final feature differentiates a line-drive type kick, where one would expect a lot of that velocity to carry through the bounce, from a high kick, where there will be less of a bias in one particular direction.","metadata":{}},{"cell_type":"markdown","source":"# 5. Optimal Punt Location\n\nGiven the models above, the goal is to find the landing coordinates for the football that maximizes the x-position after averaging over the returner's decisions.  The value of each landing position is given by the expectation value for the x-position as:\n\\begin{align}\n\\text{Value} &= x_\\text{land} - R_\\text{len} P(R) + P(NF)E(X|NF)\\\\\nE(X|NF) &= P(\\text{Touchback})(90 - x_\\text{land}) + (1-P(\\text{Touchback}))E(X|\\neg\\text{Touchback}), \\nonumber\n\\end{align}\nwhere $P(R), P(NF)$ result from the model in Section 2, $R_\\text{len}$ is the output from the model in Section 3, and $P(\\text{Touchback}), E(X|\\neg\\text{Touchback})$ are the calculated from the bouncing statistics in Section 4.","metadata":{}},{"cell_type":"code","source":"ens_df = pd.read_csv('/kaggle/input/exampleoptloc/ExampleOptimalLocation.csv')\n\nbounce_pred = ens_df.iloc[:,-7:-2].values\nmu_v = bounce_pred[:,0]\nmu_vp = bounce_pred[:,1]\n\ncos = ens_df['cos'].values\nsin = ens_df['sin'].values\n\n# Roughly approximate the expected x,y values of the bounces\nens_df['x_exp_bounce'] = mu_v*cos - mu_vp*sin\nens_df['y_exp_bounce'] = mu_v*sin + mu_v*cos\n\n# Additional statistics about the punt and bounce\nens_df['x_to_ez'] = 110 - ens_df['x_land']\nens_df['y_to_oob'] = (80/3-np.abs(80/3 - ens_df['y_land'])) * np.sign(ens_df['y_land']-80/3)\nens_df['y_rat'] = ens_df['y_to_oob']/ens_df['y_exp_bounce']\nens_df['x_exp_notb'] = ens_df['x_exp_bounce']*ens_df['y_rat']\nens_df.loc[(ens_df['y_rat']>1)|(ens_df['y_rat']<0),'x_exp_notb'] = ens_df['x_exp_bounce']\nens_df['p_tb'] = 0\nens_df['ex_ntb'] = (1-ens_df['p_tb'])*ens_df['x_exp_notb']\nens_df['Field_Val'] = ens_df['x_land'] - ens_df['R_len']*ens_df['R_avg'] + ens_df['NF_avg']*(ens_df['p_tb']*(90-ens_df['x_land']) + ens_df['ex_ntb'])\nens_df['has_calcd'] = False\n\n\n\n# Bounce model outputs statistics of the bivariate normal distributions\n# Functions to be used in calculating the probabilities and expected values\n\n# Standardize the random variables in the exponential of distributions\ndef standardize(x,mu,sigma):\n    return (x-mu)/sigma\n\n# Given a set of distribution parameters, return the bivariate normal distribution function\ndef get_bivnorm(mx,my, sx,sy, r):\n    def biv_norm(x,y):\n        x_stand = standardize(x,mx,sx)\n        y_stand = standardize(y,my,sy)\n        return 1/(2*np.pi*sx*sy*np.sqrt(1-r**2)) * np.exp(-1/(2*(1-r**2))*(x_stand**2+y_stand**2+2*r*x_stand*y_stand))\n    return biv_norm\n\n# Rotate on-field coordinates to punt-direction coordinates before acquiring PDF\ndef get_rot_pdf(mx,my, sx,sy, r, ctheta,stheta):\n    pdf = get_bivnorm(mx,my, sx,sy, r)\n    def rotated_biv_norm(x,y):\n        x_rot = x*ctheta + y*stheta\n        y_rot = y*ctheta - x*stheta\n        return pdf(x_rot,y_rot)\n    return rotated_biv_norm\n\n# Calculate the expected value of x-field coordinate\ndef get_rot_ExpX(mx,my, sx,sy, r, ctheta,stheta):\n    pdf = get_bivnorm(mx,my, sx,sy, r)\n    def rotated_expX(x,y):\n        x_rot = x*ctheta + y*stheta\n        y_rot = y*ctheta - x*stheta\n        return x * pdf(x_rot,y_rot)\n    return rotated_expX\n\n# Calculate expected value of x-field coordinate when the ball bounces out of bound\n# Assuming a straight path from the bounce location to the final location, the x-value is truncated when the ball leaves play\ndef get_rot_truncatedExpX(mx,my, sx,sy, r, ctheta,stheta, y_thres):\n    pdf = get_bivnorm(mx,my, sx,sy, r)\n    def rotated_truncexpX(x,y):\n        x_rot = x*ctheta + y*stheta\n        y_rot = y*ctheta - x*stheta\n        return x*y_thres/y * pdf(x_rot,y_rot)\n    return rotated_truncexpX\n\n# Returns y-field integration bounds depending on whether the punt is closer to the y = 0 or y = 160/3 sidelines\ndef y_oob(y_togo):\n    ranges = [[-np.inf,y_togo],[y_togo,np.inf]]\n    if y_togo>0:\n        return ranges\n    return [ranges[1],ranges[0]]\n\n# Return x-field integration bounds, for integrating both touchbacks and out of bounds\ndef get_XBounds(x_ez,y_oob):\n    slope = x_ez/y_oob\n    def tb_xbound(y):\n        x_minbound = slope*(y-y_oob)+x_ez\n        return [x_minbound, np.inf]\n\n    def oob_xbound(y):\n        x_upbound = slope*(y-y_oob)+x_ez\n        return [x_ez, x_upbound]\n    return [tb_xbound,oob_xbound]\n\n\n\n# For each punt, loop over highest field value locations and perform the full integration of expected values\n# If the Field Value remains maximal, the location is the true maximum field value\n# This works since we over-approximating the Field Value with earlier approximation\n# When the full evaluation remains maximal, or we arrive at a point which we have already calculated, we are done\n\nhighest_newVal = 0\nconverge = False\nwhile not converge:\n    best_id = ens_df.sort_values('Field_Val',ascending=False).head(1).index[0]\n\n    xl = ens_df.loc[best_id,'x_land']\n    yl = ens_df.loc[best_id,'y_land']\n\n    x_to_ez = ens_df.loc[best_id,'x_to_ez']\n    y_to_oob = ens_df.loc[best_id, 'y_to_oob']\n\n    # Only evaluate the integral if we have not explicitly integrated this point already\n    if not ens_df.loc[best_id, 'has_calcd']:\n        y_bounds = y_oob(y_to_oob)\n        x_bound_tb, x_bound_oob = get_XBounds(x_to_ez, y_to_oob)\n\n        # Get relevant functionals\n        rotated_PDF = get_rot_pdf(*bounce_pred[best_id], cos[best_id], sin[best_id])\n        rotated_ExpX = get_rot_ExpX(*bounce_pred[best_id], cos[best_id], sin[best_id])\n        rot_trunc_ExpX = get_rot_truncatedExpX(*bounce_pred[best_id], cos[best_id], sin[best_id],y_to_oob)\n\n        # Evaluate the integrals\n        # Note, it is faster to split the conditional integration bounds due to the decrease in function calls\n        ex_ip = nquad(rotated_ExpX, [[-110,x_to_ez],y_bounds[0]])[0]\n        p_tb_new = nquad(rotated_PDF, [[x_to_ez,np.inf],y_bounds[0]])[0] + nquad(rotated_PDF, [x_bound_tb,y_bounds[1]])[0]\n        ex_oob = nquad(rot_trunc_ExpX,[[-110,x_to_ez],y_bounds[1]])[0] + nquad(rot_trunc_ExpX,[x_bound_oob,y_bounds[1]])[0]\n\n        # When evaluating a point near the sideline, can update touchback probabilities in the middle of the field\n        # Unless already explicitly evaluated those locations\n        # The logic is identical to the previous step\n        # Note, this is only a useful speed-up when punting near the endzone\n        interior_locations = ens_df.loc[\n            (ens_df['x_land']==xl) \\\n            & (~ens_df['has_calcd']) \\\n            & (ens_df['y_land'].between(*sorted([80/3, ens_df.loc[best_id, 'y_land']]),inclusive='both'))].index.values\n\n        # Update the integrated values for locations closer to the middle of the field than the current location\n        ens_df.loc[interior_locations, 'p_tb'] = p_tb_new\n\n        # E(X|no TB) * P(no TB) unsure whether it should generally increase as we move towards from the middle of the field\n        # Therefore, do not update this value except when exactly calculating it\n        ens_df.loc[best_id, 'ex_ntb'] = (ex_ip + ex_oob)\n\n        # Update expected value of field position for any interior point\n        ens_df.loc[interior_locations,'Field_Val'] = ens_df.loc[interior_locations,'x_land'] \\\n                - ens_df.loc[interior_locations,'R_avg']*ens_df.loc[interior_locations,'R_len'] \\\n                + ens_df.loc[interior_locations,'NF_avg']*( \\\n                            ens_df.loc[interior_locations,'p_tb']*(90-ens_df.loc[interior_locations,'x_land']) \\\n                            + ens_df.loc[interior_locations,'ex_ntb'])\n\n    ens_df.loc[best_id, 'has_calcd'] = True\n    if best_id == ens_df.sort_values('Field_Val',ascending=False).head(1).index[0]:\n        converge = True\n        \n# This code provides the algorithm needed to calculate the optimal field location.\n# Here, it is applied to only one punt, but the same algorithm can be used to generate all of the optimal landing locations for league-average punts.\n# For the sake of time, we have included these a data set for all the 2020 punts, which we will use below\n\noptloc = pd.read_csv('/kaggle/input/fieldvalue2020/PunterOptLoc2020.csv')\noptloc['FVpastLoS'] = optloc['Field_Val'] - optloc['x_fb']\n\nfig, ax = plt.subplots(2,1,figsize=(10,8),sharex=True)\noptloc.plot.scatter(x='x_fb',y='FVpastLoS', ax=ax[0])\noptloc.plot.scatter(x='x_fb',y='p_tb', ax=ax[1])\nplt.xlabel('Line of Scrimmage (Absolute Yardline)')\nax[0].set_ylabel('Average Yards Gained')\nax[1].set_ylabel('Touchback Probability')\nplt.xlim(10,80)\nax[0].set_ylim(25,55)\nax[1].set_ylim(0,0.35)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-06T17:42:48.204667Z","iopub.execute_input":"2022-01-06T17:42:48.204999Z","iopub.status.idle":"2022-01-06T17:45:55.998052Z","shell.execute_reply.started":"2022-01-06T17:42:48.204968Z","shell.execute_reply":"2022-01-06T17:45:55.997101Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using the 2020 punt plays, the optimal landing location achievable by a league-average punter is calculated.  The average yards gained by such a punt remains relatively constant in the middle of the field, although there is a larger variance in these punts due to return teams setting up for returns.  Once the line of scrimmage is past the punt team's own 45 yardline (i.e. 55 absolute yardline), the punts are at risk of bouncing into the endzone as returners tend to let the ball bounce in hopes of a touchback.  However, the touchback probability tends to decrease during this region as the punt can be kicked higher allowing the gunners to stop the ball from entering the endzone.  The optimal strategy adds to this effect by placing the ball near the sidelines to facilitate the ball bouncing out of bounds before bouncing into the endzone.","metadata":{}},{"cell_type":"markdown","source":"# 6. Comparison of Punters in 2020\n\nEach individual punt with the punt statistics can be evaluated using the models above to find the average expected value of the punt.  In principle, the difference between the average expected value and the actual value can be used to relative values to the gunners, jammers, or returners.  Instead, comparing the average expected value of the actual punt to the optimal location for the league-average punter (labeled ΔFV) can be used to score the performance of each punter's decision making and control of the football.","metadata":{}},{"cell_type":"code","source":"players = pd.read_csv(os.path.join(dir_name,'players.csv'))\nid_to_name_dict = players.set_index('nflId')['displayName'].to_dict()\n\nreal_val = pd.read_csv('/kaggle/input/fieldvalue2020/RealLocationValues.csv')\nreal_val['Field_Val'] = real_val['x_land'] - real_val['R_avg']*real_val['R_len'] + real_val['NF_avg']*(real_val['p_tb']*(90-real_val['x_land'])+real_val['ex_ntb'])\n\n# Balls kicked directly out of bounds or into the endzone have their values set\nreal_val.loc[(real_val['y_land']<0)|(real_val['y_land']>160/3), 'Field_Val'] = real_val.loc[(real_val['y_land']<0)|(real_val['y_land']>160/3), 'x_land']\nreal_val.loc[real_val['x_land']>110,'Field_Val'] = 90\n\ncompare_rv_opt = real_val.merge(optloc[['gameId','playId','x_land','y_land','p_tb','Field_Val','kickerId']], on=['gameId','playId'], suffixes=['', '_opt'])\ncompare_rv_opt['dFV'] = compare_rv_opt['Field_Val'] - compare_rv_opt['Field_Val_opt']\n\nnum_punts = compare_rv_opt.groupby('kickerId').count()['gameId']\ntotal_diff = compare_rv_opt.groupby('kickerId')['dFV'].sum()\n\npunter_stats = pd.concat([total_diff.div(num_punts),num_punts],axis=1).reset_index().rename({0:'ΔFV Total', 'gameId':'num_punts'},axis=1)\n\npunter_vel = real_val.merge(punter[['gameId','playId','x','y']]).merge(all_pos[['gameId','playId','hangTime']]).copy()\npunter_vel['dist'] = np.sqrt((punter_vel['x']-punter_vel['x_land'])**2 + (punter_vel['y']-punter_vel['y_land'])**2)\npunter_vel['v_punt'] = calc_const(punter_vel['dist'], punter_vel['hangTime'])\npunter_stats = punter_stats.merge(punter_vel[(punter_vel['y_land']>0)&(punter_vel['y_land']<160/3)&(punter_vel['x_land']<110)].merge(optloc[['gameId','playId','kickerId']]).groupby('kickerId')['v_punt'].mean().round(2).reset_index())\npunter_stats['ΔVel'] = (punter_stats['v_punt'] - vc_avg).round(2)\n\nmidfield_punts = compare_rv_opt[compare_rv_opt['x_fb']<50].copy()\npunter_stats = punter_stats.merge(midfield_punts.groupby('kickerId')['dFV'].sum().div(midfield_punts.groupby('kickerId').count()['gameId']).reset_index().rename({0:'ΔFV Mid Field'},axis=1))\n\npinning_punts = compare_rv_opt[compare_rv_opt['x_fb']>=50].copy()\npunter_stats = punter_stats.merge(pinning_punts.groupby('kickerId')['dFV'].sum().div(pinning_punts.groupby('kickerId').count()['gameId']).reset_index().rename({0:'ΔFV Deep'},axis=1))\n\npunter_stats['Name'] = punter_stats['kickerId'].map(id_to_name_dict)\npunter_stats = punter_stats.sort_values('ΔFV Total',ascending=False).set_index('Name')\npunter_stats = punter_stats[['kickerId', 'ΔFV Total', 'ΔFV Mid Field', 'ΔFV Deep', 'ΔVel', 'v_punt','num_punts']]\n\ntop_perf = punter_stats.iloc[:,1:5].style.background_gradient(cmap='coolwarm', subset=pd.IndexSlice[:, ['ΔFV Total', 'ΔFV Mid Field', 'ΔFV Deep']])\\\n     .background_gradient(cmap='coolwarm', subset=pd.IndexSlice[:, 'ΔVel'], vmin=-1.5, vmax=1.5)\n\ndisplay(top_perf)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-06T22:07:46.733324Z","iopub.execute_input":"2022-01-06T22:07:46.733687Z","iopub.status.idle":"2022-01-06T22:07:46.850324Z","shell.execute_reply.started":"2022-01-06T22:07:46.733591Z","shell.execute_reply":"2022-01-06T22:07:46.848969Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For the 2020 data, the average difference in expected field position when compared to the optimal location for an average punter is calculated in three ranged: When all of the punts are included (ΔFV Total), punts where the original line of scrimmage is before the punt team's 40 yardline (ΔFV Mid Field), and when the line of scrimmage is past the punt team's 40 yardline (ΔFV Deep).  A better performance indicates a larger value (positive values, or least negative values).  The values of ΔFV Deep are generally much lower since optimal punts have relatively low probabilities of a touchback.  However being on the low side of this metric indicates much higher than expected rates for touchbacks.","metadata":{}},{"cell_type":"markdown","source":"### Appendix\n\nFull code and instructions to reproduce the models above can be found at: github.com/rxsims/NFLBigDataBowl2022","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-01-06T09:45:51.59891Z","iopub.execute_input":"2022-01-06T09:45:51.59923Z","iopub.status.idle":"2022-01-06T09:45:51.624789Z","shell.execute_reply.started":"2022-01-06T09:45:51.599198Z","shell.execute_reply":"2022-01-06T09:45:51.623843Z"},"_kg_hide-output":true,"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]}]}