{"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":"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(r'/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":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:42:48.278887Z","iopub.execute_input":"2022-11-09T11:42:48.279467Z","iopub.status.idle":"2022-11-09T11:42:48.320644Z","shell.execute_reply.started":"2022-11-09T11:42:48.279357Z","shell.execute_reply":"2022-11-09T11:42:48.319746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Summary\nContrary to previously published data, my analysis indicates that neither the draft percentage (staying in a position to minimise aerodynamic drag) or the path length, taking a shorter total path than the competitors, positively correlated with either the average speed or finishing place. The parameters that correlated positively both with the average speed during the race and the position at finish were related to the spatial variation of the instantaneous speed during the race, indicating that the best strategy for a racing horse is to achieve an optimal speed curve along the track (Mercier and Aftalion 2020). The results recommend that jockeys should not pursue risky racing strategies such as being close to the rails in the turns or staying too close behind other competitors to maximize draft and minimize aerodynamic drag, strategies that can increase the risk of accidents. The jockeys can in turn strive for achieving the optimal speed curve for the race and the horse. Applying XGBoost regression to the data obtained by extracting features related to the speed variation during a race, I was able to predict the position at finish with accuracy of 84% as compared to 20%, the baseline accuracy obtained by using the winning odds alone. ","metadata":{}},{"cell_type":"markdown","source":"### Main message:\nAerodynamic drafting and path efficiency do not positively correlate with winning the race or with the average speed during the race. <br>\nJockeys **shoud not** target a \"strategic race\", in the sense of minimizing aerodynamic drag or total distance ran (by running close to the rails. Jockeys should strive for a smooth race with minimum speed deviations and small drops at the turns.  \n","metadata":{}},{"cell_type":"markdown","source":"# Theoretical considerations","metadata":{}},{"cell_type":"markdown","source":"\n\n### Conclusion for this part. The speed during racing should be approximately constant at the stretches and constant at the turns, albeit with different speed values for the two race parts.  \n\nIt is not necessary to fully understand this part. If you prefer, you can skip and only remember the conclusion: the \"optimal\" race should be at constant speed at the stretch, smaller, constant speed at the turn and a comeback to the same initial speed at the stretch. This velocity profile minimizes the power used during the race. \n\nLet's start with some \"pen and paper\" estimations so that we know what to look for. The energy spent during racing is made out of three parts: external energy, correspondiong to the center of mass (CM) of the horse movement, internal energy, that is the kinetic energy of horse's body elements while they move with respect to the CM, and the energy loss due to aerodynamic drag (Davies, Spence and Wilson 2019; Pfau , Witte and M Wilson 2006). \nThe kinetic energy of the horse is $E_{kinetic} = \\frac{mv^{2}}{2}$ where $m$ is the horse's mass and $v$ is the horse's speed. Because during racing there is constant shifting from kinetic to potential energy and the internal energy is a proportion of the external energy, the internal and external energy of the galloping horse will be proportional to the kinetic energy, that is, $E_{ie} = \\alpha mv^{2}$. The power is energy spent per unit time. The power to keep the horse at a steady state of velocity $v$, by dimensional analysis is $\\frac{E_{ie}}{\\tau}$ where $\\tau$ is a system characteristic time. The most natural characteristic time is the ratio between the stride length $s$ and the average velocity, that is, $\\tau = \\frac{s}{v}$, therefore the power expenditure at each instant is $P_{ie} = \\alpha \\frac{mv^{3}}{s}$ where $\\alpha$ is a constant. Alternatively, $P_{ie} \\propto \\frac{mv^{3}}{s} $ where the $\\propto$ symbol means proportional to. The power expenditure $P_{ie}$ corresponds to the internal and external energy expenditure of galloping. This result can also be obtained via a simple model with the horse jumping at each stride for the stride length and getting an \"energy burst\" at each landing to compensate for the energy lost when impacting the soil: the total energy is $k \\cdot \\frac{mv^{2}}{2}$ because the kinetic energy and potential energy convert continuously into each other during the run. The constant $k$ takes into consideration the existence of the \"internal\" energy. At each contact with the ground, there is a fractional energy loss with a fraction $f<1$. The horse has to compensate this loss so that it continues running at the set speed. It has to \"inject\" into the race a quantity of $fk \\frac{mv^{2}}{2}$ at each contact with the ground, that is, each $\\frac{s}{v}$, where $s$ is the stride length and $v$ is the velocity. Putting them together we get for the power expenditure $fk \\frac{mv^{3}}{s}$. Now we can identify our constant $\\alpha = fk$. I know this looks like a sort of \"spherical cow\" type of physicist estimation, but keep in mind that we want a qualitative answer to the following question: should the horse vary its speed during a race? The answer is below.\n \nThe aerodynamic drag also follows a $v^{3}$ dependence (Spence et al. 2012): $P_{ad} = \\frac{1}{2}C_{d}\\rho A v^{3}$ where $C_{d}=0.9$ is an adimensional constant, $\\rho$ is the density of air, and $A$ is an area of $1m^{2}$ \n\nAdding the power due to drag and the power due to the internal and external energy, we see that the total instantaneous power expenditure $P_{t} = P_{ie}+P_{ad}$  is proportional to $v^{3}$, $P_{t} = (\\alpha \\frac{m}{s} +\\frac{1}{2}C_{d}\\rho)v^{3}$, or $P_{t} \\propto v^{3}$. The result that the total instantaneous power expenditure for a race on a straght line is proportional to the speed to the third power is new to my knowledge. Until now only the power expenditure due to the air drag has been identified as being proportional to the speed at the third power. \n\nNow that we have this proportionality relationthip we can make the following observation: to minimize the average power for a certain path with the same average velocity, the velocity along the path should be constant. The function $v^{3}$ is a convex function, therfore $P(\\frac{v_{1}+v_{2}}{2}) < \\frac{1}{2}(P(v_{1})+P\\left(v_{2})\\right)$ meaning that for the same average velocity the path with minimum power is the one where the velocity is constant. This insight can be checked numerically, see Appendix. Obviously the velocity cannot be constant in a horserace because all horses start from rest, but our analysis shows that after the initial acceleration phase, for a straight course, the optimal strategy is to run at constant speed. In (Spence et al. 2012) there is a mention that \"for middle distance events (greater than 2 minutes) constant pace may be optimal\" but there is no argument as of why. Our data show that even in short races, the individual horses cannot keep the pace for the whole race. In (Mercier and Aftalion 2020 ) where each horse was analysed alone, the hypothesis was that there is an \"optimal\" speed vs. distance dependence for each horse. \n\nAt the turn we have the centripetal force $F_{c}=\\frac{mv^{2}}{R}$, where $R$ is the radius of the turn. We cannot simply use the well known function for the power $P=F\\cdot v$ because the vectors of the centripetal force and the velocities are perpendicular here. Still, there is an extra metabolic consumption, think of isostatic effort (such as keeping a heavy load) that requires metabolic expenditure even though the mechanical work is zero (due to zero displacement). As a dimensional analysis first indication, the power due to the centripetal force is $P_{c} \\propto v^{3}$ or $P_{c} = \\beta v^{3}$, therefore the total instantaneous power in the turn is $P_{t\\_turn} (\\alpha \\frac{m}{s} +\\frac{1}{2}C_{d}\\rho + \\beta)v^{3}$ or $P_{t\\_turn}  \\propto v^{3}$ just with a larger proportionality constant than the case of the straght line. \n\nOur simple analysis indicates that the speed during a race should be constant at the stretch, then jump to a different (lower) speed at the turn and then constant during the last stretch. Our estimation and data also show that the goal of this sort of \"almost constant\" pace strategy should be pursued for short races also. \n ","metadata":{}},{"cell_type":"markdown","source":"### Percentage effects of air drag and path efficiency are expected to be small\n\nExisting literature puts the contribution of air drag to the total energy expenditure during the race at 17%. (Davies, Spence and Wilson 2019). If the air drag is reduced by a factor of say, 20%, due to \"strategic riding\", the energy saved would be around 3.4%. For the estimation of the air drag the horse and jockey are assumed a bluff body of one meter squared. Because the horse's shape is aerodynamical and it is not really a bluff body, chances are that the value obtained above for the drag is an upper limit and the actual value is lower. Experiments in wind tunnels are needed to measure exactly the air drag. \nPath efficiency can be another way to save energy during the race. At the turns the horses farther away from the rails run longer distances and it is interesting to see whether being closer to the rails and minimizing the distance ran brings about any benefit in terms of average velocity and chance of winning. \n\nNotice that trying to run close to the rails or running in such a way to optimize air drag is risky: close to the rails the track quality is lower, the horse may get boxed in, and running close to another competitor or in a compact group increases the risk of accidents. Therefore it is of great interest, for the health and safety understanding within the sport, to know whether pursuing such \"strategic riding\" as in minimizing running distance of minimizing air drag are worth it. ","metadata":{}},{"cell_type":"markdown","source":"### What hypothesis will we check? \nIn the literature available there are two lines of thought: look at the influence of the other horses in the race, especially at minimizing air drag by sitting in a \"drafting position\" behind another competitor (Spence et al. 2012).\nAnother view is to look at each horse run independently and identify some sort of \"ideal\" velocity profile along the track (Mercier and Aftalion 2020). \n\nThe question is: what should a jockey do during a race? Try to minimize air drag and distance ran, or try to optimize the speed profile during the race? The hypothesis can be for example \"the speed profile during the race is more important than the drafting\". It is reasonable, can be experimentally disporven (negated), therefor it is a scientific hypothesis. \n\nTo answer this question, the plan is to measure the effects of the drafting, the path efficiency and the speed profiles. According to the definitions in Physics college classes, we will use from now on velocity for the velocity vector and speed for the modulus of the velocity, a scalar.   ","metadata":{}},{"cell_type":"markdown","source":"# Materials and Methods","metadata":{}},{"cell_type":"markdown","source":"We begin by transforming the latitude and longitude coordinates into $x$ and $y$ coordinates. The $x$ and $y$ coordinates will then be rotated into $x_{r}$ and $y_{r}$ and and translated so that the chute is at negative $x_{r}$ values, the chute and the first stretch are paralel to the $x_{r}$ axis and the center of the leftmost turn (the one above the chute) is at position $ x_{r}= 0$ and $y_{r}=0$. With these coordinates axes it will be straightforward to find the actual projected position (the position along the track) for each point in the plane, that is, for each racing horse. The actual projected position is the position that counts in terms of advancing along the track during the race. \n","metadata":{}},{"cell_type":"code","source":"!pip install pymap3d\n!pip install mpl-scatter-density\nimport pymap3d as pm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport math\nfrom numpy import linalg as LA\nfrom numpy.linalg import inv\nimport torch\nfrom scipy import  polyfit\nfrom scipy import optimize\nfrom numpy import *\nfrom scipy.signal import savgol_filter\nfrom scipy import stats\nfrom mpl_toolkits.mplot3d import Axes3D\nimport plotly.express as px\nfrom tqdm import tqdm\nplt.style.use('ggplot')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:42:48.373848Z","iopub.execute_input":"2022-11-09T11:42:48.374803Z","iopub.status.idle":"2022-11-09T11:43:20.290868Z","shell.execute_reply.started":"2022-11-09T11:42:48.374763Z","shell.execute_reply":"2022-11-09T11:43:20.289529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\nimport xgboost as xgb\nfrom scipy import stats\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import accuracy_score\nimport shap\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import GridSearchCV\nimport sklearn\nimport warnings # Supress warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:20.293670Z","iopub.execute_input":"2022-11-09T11:43:20.294139Z","iopub.status.idle":"2022-11-09T11:43:22.591406Z","shell.execute_reply.started":"2022-11-09T11:43:20.294096Z","shell.execute_reply":"2022-11-09T11:43:22.589333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BASE_DIR=r'/kaggle/input'","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:22.593183Z","iopub.execute_input":"2022-11-09T11:43:22.593564Z","iopub.status.idle":"2022-11-09T11:43:22.599269Z","shell.execute_reply.started":"2022-11-09T11:43:22.593530Z","shell.execute_reply":"2022-11-09T11:43:22.597811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### We put all data files in a dict and we decide on what type of races to focus","metadata":{}},{"cell_type":"code","source":"files_dict = {}\nfor dirname, _, filenames in os.walk(r'/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n        files_dict[filename] = pd.read_csv(os.path.join(dirname, filename))\n\nfiles_list=[]\nfor key in files_dict.keys():\n    print(key)\n    files_list.append(key)\nprint(files_list)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:22.604545Z","iopub.execute_input":"2022-11-09T11:43:22.605047Z","iopub.status.idle":"2022-11-09T11:43:50.581928Z","shell.execute_reply.started":"2022-11-09T11:43:22.605008Z","shell.execute_reply":"2022-11-09T11:43:50.580745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race = files_dict['nyra_race_table.csv']\ndf_race['furlongs'] = df_race['distance_id']/100","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:50.594697Z","iopub.execute_input":"2022-11-09T11:43:50.595642Z","iopub.status.idle":"2022-11-09T11:43:50.608888Z","shell.execute_reply.started":"2022-11-09T11:43:50.595580Z","shell.execute_reply":"2022-11-09T11:43:50.607699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### I will look at short races, that is races 7 furlongs and less in length. \nThis choice fits with previously published results, on a much larger dataset where they looked at short races, medium range races and long races (Spence et al. 2012). In our dataset the short races are the most numerous,  more than half the total number of races are short races and I will focus on them for now. The number of long races is simply too small compared to the number of short races and to the number of previously analysed races to get a good indication of what factors influence outcome. ","metadata":{}},{"cell_type":"code","source":"df_race.hist('furlongs', bins = 20)\nplt.title('The majority of the races are shorter than 7 furlongs.')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-09T11:43:50.611710Z","iopub.execute_input":"2022-11-09T11:43:50.612577Z","iopub.status.idle":"2022-11-09T11:43:50.934822Z","shell.execute_reply.started":"2022-11-09T11:43:50.612522Z","shell.execute_reply":"2022-11-09T11:43:50.933580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The graph above also shows very few \"long races\" with distance larger than 2890m or 14 furlongs, and in the \"medium races\" there is a peak at the 8-9 furlongs distance with a tail up to 12 furlongs. The more numerous group of races is the \"short races\" and this is where we are going to look further. ","metadata":{}},{"cell_type":"markdown","source":"Following (Spence et al. 2012) we are going to analyse sprints in our work, that is races 7 furlongs and less. ","metadata":{}},{"cell_type":"code","source":"long_races = df_race.query('furlongs > 7.05')\nnumber = len(long_races)\npercentage = len(long_races)/len(df_race)\nprint(f'The percentage of long races is {percentage: .1%}.' )\nprint(f'The number of long races is {number}.' ) \n                       ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:50.936851Z","iopub.execute_input":"2022-11-09T11:43:50.937727Z","iopub.status.idle":"2022-11-09T11:43:50.952144Z","shell.execute_reply.started":"2022-11-09T11:43:50.937661Z","shell.execute_reply":"2022-11-09T11:43:50.950721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"short_races = df_race.query('furlongs< 7.05')\npercentage = len(short_races)/len(df_race)\nnumber = len(short_races)\nprint(f'The percentage of short races is {percentage: .1%}.' )\nprint(f'The number of short races is {number}.' )                     ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:50.953520Z","iopub.execute_input":"2022-11-09T11:43:50.954397Z","iopub.status.idle":"2022-11-09T11:43:50.966785Z","shell.execute_reply.started":"2022-11-09T11:43:50.954343Z","shell.execute_reply":"2022-11-09T11:43:50.965163Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### I have a fit for 8 of the racetracks on the three venues (Aqueduct Inner turf is left out). \nEach track on the racetrack is composed of two line segments and two semicircles. We can fit them to points collected by hand. You will see said segments and said semicircles when I plot them after the fit.  I have collected by hand using google maps a collection of points for each segment and each semicircle for each track except Aqueduct inner turf where the rails were not visible in the satellite data. The segments and semicircle points are collected in the attached csv files with self explanatory names: aqu_d_c1 means the points from the semicircle c1 on Aqueduct dirt \"course_type\". \n\nEach racetrack has three tracks, the D, the turf T and the inner turf I. The names for the Aqueduct are different, D, O and I, and we will take that aspect into consideration. In the code below I load the points collected by hand so that I can use them to get back the full fit for the track. This way I will have for each venue and course type analyzed an exact indication where the rails are. ","metadata":{}},{"cell_type":"markdown","source":"### Here I load the complete dataset, inlcuding the points for the fits, in a dataframe called df_complete","metadata":{}},{"cell_type":"code","source":"df_complete = files_dict['nyra_2019_complete.csv'].copy()\ndf_complete = df_complete.rename(columns={\n    'AQU': 'track_id', \n    '2019-01-01':'race_date', \n    '9': 'race_number',\n    '6  ':'program_number',\n    '72': 'trakus_index',\n    '40.6729017197787': 'latitude',\n    '-73.8276065972899': 'longitude',\n    '600': 'distance_id',\n    'D': \"course_type\",\n    'GD ': 'track_condition',\n    '48': 'run_up_distance',\n    'CLM': 'race_type',\n    '25000.00': 'purse',\n    '00420': 'post_time',\n    '120': 'weight_carried',\n    'Andre Shivnarine Worrie': 'jockey',\n    '2090': 'odds', \n    '8': 'position_at_finish'}\n    )","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:50.968520Z","iopub.execute_input":"2022-11-09T11:43:50.969230Z","iopub.status.idle":"2022-11-09T11:43:52.432478Z","shell.execute_reply.started":"2022-11-09T11:43:50.969188Z","shell.execute_reply":"2022-11-09T11:43:52.431152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# here I just eliminate the spaces from the program numbers strings\ndf_complete['program_number'] = df_complete['program_number'].astype('str')\ndf_complete['program_number'] = df_complete['program_number'].str.strip()\ndf_complete.program_number.unique()","metadata":{"execution":{"iopub.status.busy":"2022-11-09T11:43:52.436407Z","iopub.execute_input":"2022-11-09T11:43:52.436827Z","iopub.status.idle":"2022-11-09T11:43:55.229391Z","shell.execute_reply.started":"2022-11-09T11:43:52.436789Z","shell.execute_reply":"2022-11-09T11:43:55.228173Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_complete.course_type.unique()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:55.230872Z","iopub.execute_input":"2022-11-09T11:43:55.231247Z","iopub.status.idle":"2022-11-09T11:43:55.561998Z","shell.execute_reply.started":"2022-11-09T11:43:55.231212Z","shell.execute_reply":"2022-11-09T11:43:55.560714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here I load the points taken by hand, points used to fit every track for the racetracks and add them to the dataframe. ","metadata":{}},{"cell_type":"code","source":"all_points = pd.read_csv(r'/kaggle/input/points-for-fitting-tracks/all_points.csv')\n# here I add the points that I will use to fit the racetracks \ndf_complete = pd.concat([df_complete, all_points])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:55.563905Z","iopub.execute_input":"2022-11-09T11:43:55.564285Z","iopub.status.idle":"2022-11-09T11:43:57.030792Z","shell.execute_reply.started":"2022-11-09T11:43:55.564253Z","shell.execute_reply":"2022-11-09T11:43:57.029663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I also add a distance column and a winning probability column to the data. ","metadata":{}},{"cell_type":"code","source":"df_complete['distance'] = (df_complete['distance_id']/100)*201.168 + df_complete['run_up_distance']*0.3048\ndf_complete['odds'] = df_complete['odds']/100\ndf_complete['win_prob'] =  1/(df_complete['odds'] + 1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:57.032408Z","iopub.execute_input":"2022-11-09T11:43:57.032813Z","iopub.status.idle":"2022-11-09T11:43:57.204162Z","shell.execute_reply.started":"2022-11-09T11:43:57.032775Z","shell.execute_reply":"2022-11-09T11:43:57.203152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Following Joe Appelbaum's advice in the discussions I will look only at course_type D, O, I, T and leave M aside. \nAs seen from the pie plot below, we still keep most of the races. ","metadata":{}},{"cell_type":"code","source":"plt.title('Course types pie')\nplt.pie(df_complete['course_type'].value_counts(), labels=df_complete['course_type'].unique())\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:57.205689Z","iopub.execute_input":"2022-11-09T11:43:57.206427Z","iopub.status.idle":"2022-11-09T11:43:57.882481Z","shell.execute_reply.started":"2022-11-09T11:43:57.206386Z","shell.execute_reply":"2022-11-09T11:43:57.880837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### How I got the positions of the horses and tracks so that the analysis is straightforward\nThe following code transforms the latitude and longitude into x and y coordinates and then rotates and translates said (x,y) coordinates into (xr, yr) coordinates such that the leftmost semicircle of each track (D, O, I, or T) has its center at the (0,0) position. When I say \"track\" I mean the running surface for a certain course_type and certain venue. For example the D track of the AQU venue. The steps are as follows:\n> -   for a certain venue, transform all the (latitude,longitude) coordinates into (x,y,z) coordinates using the pymap3d python module. I do not worry about the  few 20m or so of \"altitude\" for the venues, the error is small when compared with the radius of the Earth and other errors in the data. \n> -  for a certain venue, make a coordinates transform that will render all points in the (x,y) plane. This step is carried out by calcuating the moments of inertia for the object created by all the position points and getting the transformation matrix from the eigenvectors of the matrix containing the moments of inertia. \nAt the end of this step, for each venue the position points are in an (x,y) plane. Now let's look at each course type in each venue. \n> - for a course type, fit the lines corresponding to the two stretches \n> - for the same course type,  rotate the axes of coordinates so that the two stretches are horizontal\n> - for the same course type, translate the axes of coordinates so that the center of the leftmost semicircle is at (0,0).\n> - for the same course type find the radius of the turn by fitting the points collected at the turn\n>  These coordinates are designates with xr, yr. \nRepeat these steps for each course types. \n","metadata":{}},{"cell_type":"markdown","source":"We defined a number of functions to do the steps described above. You do not have to understand or follow all the steps, you can look at the graph outputs in the notebook. ","metadata":{}},{"cell_type":"code","source":"def add_xy(df_venue):\n    '''\n    This function adds x and y values to the dataframe for AQU and BEL venues (at a time). For SAR there is another \n    function defined. after transforming the latitude and longitude into the ecef coordinates we obtain a set of points in x',y',z'. \n    To transform these x',y',z' axes into an xy plane a straightforward way to do it is to calculate the moments of inertia, \n    the eigenmatrices corresponding to these moment of inertia, and then use the inverses of these matrices to transform the \n    x',y',z' coordinate axes. The \"object\" whose moments of inertia we calculate is composed of all the tracks that took place on \n    the particular venue that we are finding the x and y for. \n    '''\n    \n    \n    df_venue.reset_index(inplace=True)\n    df_venue=df_venue.drop(columns=['index'])\n    ecef_T = np.array(pm.geodetic2ecef(df_venue[\"latitude\"].values, df_venue[\"longitude\"].values, np.array([2]*len(df_venue))))\n    # we find the center of mass\n    xcm=ecef_T[0].mean()\n    ycm=ecef_T[1].mean()\n    zcm=ecef_T[2].mean()\n    #these are the coordinates with respect with the center of mass\n    ecef_cm = ecef_T - np.array([[xcm]*len(ecef_T[0]), [ycm]*len(ecef_T[0]),[zcm]*len(ecef_T[0])])\n    ecef_cm_t = ecef_cm.T\n    Ixx=np.dot(ecef_cm[1],ecef_cm[1])+np.dot(ecef_cm[2],ecef_cm[2])\n    Iyy=np.dot(ecef_cm[0],ecef_cm[0])+np.dot(ecef_cm[2],ecef_cm[2])\n    Izz=np.dot(ecef_cm[0],ecef_cm[0])+np.dot(ecef_cm[1],ecef_cm[1])\n    Ixy=-np.dot(ecef_cm[0],ecef_cm[1])\n    Ixz=-np.dot(ecef_cm[0],ecef_cm[2])\n    Iyz=-np.dot(ecef_cm[1],ecef_cm[2])\n    I_mat=np.matrix([[Ixx,Ixy,Ixz],[Ixy,Iyy,Iyz],[Ixz,Iyz,Izz]])\n    w, v = LA.eig(I_mat)\n    # w are the eigenvalues, v is the matrix that has the eigenvectors as columns. it is a 3x3 matrix\n    # tr_mat is the transformation matrix, the inverse of the inverse transposed\n    tr_mat = inv(inv(v).T)\n    # in transformed the coordinates into tensors so that I can apply batch matrix multiplication to all dimension 3 \n    # x',y',z' coordinate vectors\n    x_data = torch.tensor(ecef_cm_t)\n    tr_tensor=torch.tensor(tr_mat)\n    tr_tensor_b = tr_tensor.repeat(len(ecef_cm_t),1,1)\n    new_coord=torch.bmm(tr_tensor_b, x_data.unsqueeze(dim=2))\n    coord_np=new_coord.squeeze(axis=2).numpy()\n    coord_np_t=coord_np.T\n    # when plotting them I saw wich is which. This is why for SAR another function was needed because I had to take\n    # other elements of the coord array (1,0) for SAR instead of (2,0) for BEL and AQU.\n    df_venue['x'] = coord_np_t[2]\n    df_venue['y'] = coord_np_t[0]\n    #df_venue['z'] = coord_np_t[1]\n    return df_venue","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:57.885386Z","iopub.execute_input":"2022-11-09T11:43:57.886458Z","iopub.status.idle":"2022-11-09T11:43:57.918280Z","shell.execute_reply.started":"2022-11-09T11:43:57.886394Z","shell.execute_reply":"2022-11-09T11:43:57.916643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_xy_sar(df_venue):\n    '''\n    The same as add_xy_df function above but for for saratoga. the only thing that is changed is the selection at the end. \n    '''\n    #df_venue=df_complete.query('track_id == @venue_name')\n    df_venue.reset_index(inplace=True)\n    df_venue=df_venue.drop(columns=['index'])\n    ecef_T = np.array(pm.geodetic2ecef(df_venue[\"latitude\"].values, df_venue[\"longitude\"].values, np.array([2]*len(df_venue))))\n    xcm=ecef_T[0].mean()\n    ycm=ecef_T[1].mean()\n    zcm=ecef_T[2].mean()\n    ecef_cm = ecef_T - np.array([[xcm]*len(ecef_T[0]), [ycm]*len(ecef_T[0]),[zcm]*len(ecef_T[0])])\n    ecef_cm_t = ecef_cm.T\n    Ixx=np.dot(ecef_cm[1],ecef_cm[1])+np.dot(ecef_cm[2],ecef_cm[2])\n    Iyy=np.dot(ecef_cm[0],ecef_cm[0])+np.dot(ecef_cm[2],ecef_cm[2])\n    Izz=np.dot(ecef_cm[0],ecef_cm[0])+np.dot(ecef_cm[1],ecef_cm[1])\n    Ixy=-np.dot(ecef_cm[0],ecef_cm[1])\n    Ixz=-np.dot(ecef_cm[0],ecef_cm[2])\n    Iyz=-np.dot(ecef_cm[1],ecef_cm[2])\n    I_mat=np.matrix([[Ixx,Ixy,Ixz],[Ixy,Iyy,Iyz],[Ixz,Iyz,Izz]])\n    w, v = LA.eig(I_mat)\n    # w are the eigenvalues, v is the matrix that has the eigenvectors as columns. it is a 3x3 matrix\n    tr_mat = inv(inv(v).T)\n    x_data = torch.tensor(ecef_cm_t)\n    tr_tensor=torch.tensor(tr_mat)\n    tr_tensor_b = tr_tensor.repeat(len(ecef_cm_t),1,1)\n    new_coord=torch.bmm(tr_tensor_b, x_data.unsqueeze(dim=2))\n    coord_np=new_coord.squeeze(axis=2).numpy()\n    coord_np_t=coord_np.T\n    df_venue['x'] = coord_np_t[1]\n    df_venue['y'] = coord_np_t[0]\n    new_venue = df_venue.copy()\n    return new_venue","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:57.920412Z","iopub.execute_input":"2022-11-09T11:43:57.922553Z","iopub.status.idle":"2022-11-09T11:43:57.953790Z","shell.execute_reply.started":"2022-11-09T11:43:57.922478Z","shell.execute_reply":"2022-11-09T11:43:57.951694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# this function finds the center and radius of a circle given points on its circumference. We use it \n# to fit the racetracks. \ndef fit_circle(d_x, d_y):\n    # this is a function that fits a number of points in the x y plane to a circle. I found most of the code on stackoverflow.\n    x = d_x\n    y = d_y\n    \n    def calc_R(c):\n        \"\"\" calculate the distance of each 2D points from the center c=(xc, yc) \"\"\"\n        return sqrt((x-c[0])**2 + (y-c[1])**2)\n\n    def calc_ecart(c):\n        \"\"\" calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) \"\"\"\n        Ri = calc_R(c)\n        return Ri - Ri.mean()\n    \n    x_m = x.mean()\n    y_m = y.mean()\n        \n    # calculation of the reduced coordinates\n    u = x - x_m\n    v = y - y_m\n    center_estimate = x_m, y_m\n    center_2, ier = optimize.leastsq(calc_ecart, center_estimate)\n\n    xc_2, yc_2 = center_2\n    Ri_2       = calc_R(center_2)\n    R_2        = Ri_2.mean()\n    residu_2   = sum((Ri_2 - R_2)**2)\n    return R_2, xc_2, yc_2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:57.960191Z","iopub.execute_input":"2022-11-09T11:43:57.961334Z","iopub.status.idle":"2022-11-09T11:43:57.972635Z","shell.execute_reply.started":"2022-11-09T11:43:57.961255Z","shell.execute_reply":"2022-11-09T11:43:57.971714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def rotate_shift(df_venue, course_typ):\n# This function rotates and translates the x y coordinates for each racetrack so that the (x,y)=(0,0)\n# point is at the center of the leftmost circle and the positive x axis is to the right, the \n# positive y axis is up. \n    xl1 = df_venue.query('course_type == @course_typ & race_type == \"l1\"')['x']\n    yl1 = df_venue.query('course_type == @course_typ & race_type == \"l1\"')['y']\n    xl2 = df_venue.query('course_type == @course_typ & race_type == \"l2\"')['x']\n    yl2 = df_venue.query('course_type == @course_typ & race_type == \"l2\"')['y']\n \n    (a2,b2) = np.polyfit(xl2, yl2, 1)\n    (a1, b1) =np.polyfit(xl1, yl1, 1)\n    alpha = -math.atan((a1+a2)/2)\n    # i rotate the coordinates axes so that the stretches defined by the lines l1 and l2 are parallel to the new x coordinate\n    # xr\n    df_venue['xr'] = df_venue['x']*math.cos(alpha) - df_venue['y']*math.sin(alpha)\n    df_venue['yr'] = df_venue['x']*math.sin(alpha) + df_venue['y']*math.cos(alpha)\n    \n    venue_fits = df_venue.query('race_type == [\"c1\", \"c2\",\"l1\",\"l2\"] & course_type == \"T\"')\n    #df_venue.plot(x='xr', y='yr', kind='scatter')\n    # the points of circle1\n    xc1 = df_venue.query('course_type == @course_typ & race_type == \"c1\"')['xr']\n    yc1 = df_venue.query('course_type == @course_typ & race_type == \"c1\"')['yr']\n    # the points of circle2\n    xc2 = df_venue.query('course_type == @course_typ & race_type == \"c2\"')['xr']\n    yc2 = df_venue.query('course_type == @course_typ & race_type == \"c2\"')['yr']\n\n    R1, xc1_fit, yc1_fit = fit_circle(xc1, yc1)\n    R2, xc2_fit, yc2_fit = fit_circle(xc2, yc2)\n    \n    #here  I shift the rotated axes so that the circle 1 is at (0,0)\n    df_venue['xr'] = df_venue['xr'] - xc1_fit\n    df_venue['yr'] = df_venue['yr'] - yc1_fit\n    \n    xc2 = xc2 - xc1_fit\n    xc1 = xc1 - xc1_fit\n    yc2 = yc2 - yc1_fit\n    yc1 = yc1 - yc1_fit\n    \n    theta_fit2 = linspace(-pi/2, pi/2, 180)\n    theta_fit1 = linspace(pi/2, 3*pi/2, 180)\n\n    x_fit1 =  R1*cos(theta_fit1)\n    y_fit1 =  R1*sin(theta_fit1)\n    x_fit2 = xc2_fit - xc1_fit + R2*cos(theta_fit2)\n    y_fit2 = yc2_fit - yc2_fit + R2*sin(theta_fit2)\n    \n    colors_dict = {'D': 'pink', 'T': 'thistle', 'I': 'lavender', 'O': 'plum'}\n\n    course = course_typ\n    colour = colors_dict[course_typ]\n    df_now = df_venue.query('course_type == @course')\n    x = df_now['xr'].to_numpy()\n    y = df_now['yr'].to_numpy()\n    line12x = linspace(0, xc2_fit-xc1_fit, 100)\n    line1y = [-(R1+R2)/2]*100\n    line2y = [(R1+R2)/2]*100\n    \n    plt.rcParams[\"figure.figsize\"] = (24,24)\n    plt.scatter(x, y , color=colour, label='course type: '+str(course), s=6)\n    plt.plot(x_fit1, y_fit1, 'yellow' ,  lw=6, label='circle C1')\n    plt.plot(x_fit2, y_fit2, 'r' ,  lw=6, label='circle C2')\n    plt.plot(line12x, line1y, 'blue' ,  lw=6, label='line l1')\n    plt.plot(line12x, line2y, 'green' ,  lw=6, label='line l2')\n    plt.axvline(0.0, color = 'k')\n    plt.axhline(0,0, color = 'k')\n    plt.xlabel('xr (meters)')\n    plt.ylabel('yr (meters)')\n    \n    plt.legend(loc = 'upper left')\n    plt.title('course types and course fits at ' + df_venue.track_id.unique())\n    \n    ax = plt.gca()\n    ax.yaxis.grid(True)\n    ax.xaxis.grid(True)\n    plt.gca().set_aspect('equal', adjustable='box')\n\n    # it returns the dataframe with the rotated coordinates, the distance between the circles origins and the radius\n    # of the circles\n    # -xc1_fit is the x translation vector, along the pp axis\n    return df_venue, xc2_fit - xc1_fit, (R1+R2)/2, -xc1_fit","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:57.974053Z","iopub.execute_input":"2022-11-09T11:43:57.974928Z","iopub.status.idle":"2022-11-09T11:43:57.995701Z","shell.execute_reply.started":"2022-11-09T11:43:57.974893Z","shell.execute_reply":"2022-11-09T11:43:57.994513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The adding of xr, yr coordinates, the good rotated coordinates have to be done for each venue at a time.","metadata":{}},{"cell_type":"code","source":"aqu = df_complete.query('track_id == \"AQU\" ')\naqu = add_xy(aqu)\n\nbel = df_complete.query('track_id == \"BEL\" ')\nbel = add_xy(bel)\n\nsar = df_complete.query('track_id == \"SAR\" ')\nsar = add_xy_sar(sar)\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:43:57.997020Z","iopub.execute_input":"2022-11-09T11:43:57.997369Z","iopub.status.idle":"2022-11-09T11:44:04.828468Z","shell.execute_reply.started":"2022-11-09T11:43:57.997337Z","shell.execute_reply":"2022-11-09T11:44:04.826920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In a dictionary called **dict_t** I put all the dataframes of the rotated and translated coordinates together with info on fits:  radius of the turn, length of the stretch and translation distance.","metadata":{}},{"cell_type":"code","source":"dict_t={}","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:04.830075Z","iopub.execute_input":"2022-11-09T11:44:04.830543Z","iopub.status.idle":"2022-11-09T11:44:04.837013Z","shell.execute_reply.started":"2022-11-09T11:44:04.830493Z","shell.execute_reply":"2022-11-09T11:44:04.835553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here I add xr and yr (from rotated) coordinates to each point in the AQU, BEl and SAR venues. If I plot it on x and y you will get the shape of the race track \nand the coordinates zeros x=0 and y=0 close to the center of mass. ","metadata":{}},{"cell_type":"code","source":"for j in ['D', 'T', 'I']:\n    df = bel.query('course_type == @j').reset_index(drop=True)\n    dict_t['bel'+'_'+j] = rotate_shift(df, j)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:04.838471Z","iopub.execute_input":"2022-11-09T11:44:04.838955Z","iopub.status.idle":"2022-11-09T11:44:10.615532Z","shell.execute_reply.started":"2022-11-09T11:44:04.838901Z","shell.execute_reply":"2022-11-09T11:44:10.614306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for j in ['D', 'T', 'I']:\n    df = sar.query('course_type == @j').reset_index(drop=True)\n    dict_t['sar'+'_'+j] = rotate_shift(df, j)","metadata":{"_kg_hide-input":true,"_kg_hide-output":false,"execution":{"iopub.status.busy":"2022-11-09T11:44:10.617489Z","iopub.execute_input":"2022-11-09T11:44:10.618321Z","iopub.status.idle":"2022-11-09T11:44:14.243993Z","shell.execute_reply.started":"2022-11-09T11:44:10.618270Z","shell.execute_reply":"2022-11-09T11:44:14.242778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for j in ['D', 'O']:\n    df = aqu.query('course_type == @j').reset_index(drop=True)\n    dict_t['aqu'+'_'+j] = rotate_shift(df, j)","metadata":{"_kg_hide-input":true,"_kg_hide-output":false,"execution":{"iopub.status.busy":"2022-11-09T11:44:14.245280Z","iopub.execute_input":"2022-11-09T11:44:14.245653Z","iopub.status.idle":"2022-11-09T11:44:20.141202Z","shell.execute_reply.started":"2022-11-09T11:44:14.245619Z","shell.execute_reply":"2022-11-09T11:44:20.139775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Here are the keys of the said dict', dict_t.keys())","metadata":{"execution":{"iopub.status.busy":"2022-11-09T11:44:20.142819Z","iopub.execute_input":"2022-11-09T11:44:20.143597Z","iopub.status.idle":"2022-11-09T11:44:20.150695Z","shell.execute_reply.started":"2022-11-09T11:44:20.143549Z","shell.execute_reply":"2022-11-09T11:44:20.149481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('the stretch, turn radius and translation distance for BEL turf in meters: ', dict_t['bel_T'][1:])","metadata":{"execution":{"iopub.status.busy":"2022-11-09T11:44:20.152565Z","iopub.execute_input":"2022-11-09T11:44:20.153357Z","iopub.status.idle":"2022-11-09T11:44:20.163079Z","shell.execute_reply.started":"2022-11-09T11:44:20.153320Z","shell.execute_reply":"2022-11-09T11:44:20.161814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"### The track dimensions obtained from the fits are similar to the ones indicated in the venue descriptions. \nBelow I verify that the track dimensions are what we expect from the available venue description. This analysis is not needed for my report, I just put it here in case someone wants to make sure the track dimensions I get are OK. I hid the cells and their output.  \n\nThe dimensions obtained are in good agreement with the dimensions published for the racetracks. \nhttps://www.drf.com/sites/default/files/US%20Track%20Dimensions%20List%20-%20May%202015.pdf\nFor example, for AQU Dirt, the dimensions for the stretch are given as: 1485.5 feet which is 452.7804m and we got for the stretch length 453.10m, a deviation of 0.07% or about one foot in 1485 feet. \nThe radius is given as 470 feet, that is 143.256m and we got for the two radii (they should be the same value) 143.24m and 143.095m. Again, very good agreement and we are confident that the model of the track that we got is good for our purpose of analysing the races. \nWith these axes, the \"projected position\" (pp) , that is the position that counts with respect to winning and being ahead, is easy to find. For example, for positions on the chute and the first stretch, pp = x. Do not get confused that the pp can be negative for x<0, one can have a start at a negative pp, the key is to keep correct track of the variation of pp for all races and all horses. For the turns we use the definition of the angles in radians and the pp is $xc_{2} + l\\frac{\\alpha}{\\pi}$ where $l$ is the length of turn as taken from the available data and $\\alpha$ is the angle swept measured in radians, where $xc_{2}$ is the length of the stretch, equal to the x coordinate of the second semicircle.","metadata":{"_kg_hide-input":false}},{"cell_type":"markdown","source":"BELMONT\n> Main Track: 1 1/2 Miles Around\n> 1097' Head of Stretch to Finish Line\n> 843' Finish Line to First Turn\n> 643' Turn Radius\n> 2020' Length of Turn\nA total of 2414m\nThe calculated values, extracted from my fit is below: ","metadata":{"_kg_hide-input":true}},{"cell_type":"code","source":"cal_bel = 2*dict_t['bel_D'][1]+2*dict_t['bel_D'][2]*math.pi\nknown_bel = 2414\npercent = (known_bel-cal_bel)/known_bel\nprint(f'The calculated length of the Belmont D track is {cal_bel:.2f} and the known length of the \\\n      Belmont D track is {known_bel:.2f}, a {percent: .1%} deviation')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:20.165040Z","iopub.execute_input":"2022-11-09T11:44:20.165763Z","iopub.status.idle":"2022-11-09T11:44:20.175855Z","shell.execute_reply.started":"2022-11-09T11:44:20.165725Z","shell.execute_reply":"2022-11-09T11:44:20.174492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The function below adds the \"projected\" coordinates to each track, denoted 'pp'. \nThese are the coordinates that count for being ahead. For example, if a horse is sweeping an angle $\\alpha$ while at the turn, the projected position displacement is $\\Delta pp = r\\alpha$ where $r$ is the radius of the track, while the actual position displacement is $\\Delta pa = R\\alpha$ where $R$ is the radius at which the horse is actually running, with $R > r$. At the stretch, if the real displacement is $\\Delta l$ the projected displacement is $\\Delta p = \\Delta l cos(\\alpha)$ , where $\\alpha$ is the angle the displacement vector $\\Delta \\vec{l}$ makes with the direction of the track. The projected position corresponds to the displacement that counts for being ahead and it also corresponds to the given track sizes. The zeros of the projected positons are at the beginning of the first stretch, therefore the **pp** coordinates at the chute are negative. That's OK because we look at $\\Delta pp$ all the time. ","metadata":{}},{"cell_type":"code","source":"def pp_sprint(df, stretch, r):\n    # this function adds the projected position, that is the position along the track\n    #this function works only with sprints, that is less than 7 furlongs    \n    trackus_threshold = 100\n    if df['distance_id'].max() < 710:\n        df = df.reset_index().drop(columns='index')\n        chute = df.query('xr < 0  & trakus_index < @trackus_threshold').reset_index().drop(columns = 'index')\n        stretch2 = df.query('xr > 0 & xr < @stretch & yr < 0 ').reset_index().drop(columns = 'index')\n        circle3 = df.query('xr > @stretch ').reset_index().drop(columns = 'index')\n        stretch4 = df.query('xr > 0 & xr < @stretch & yr > 0 ').reset_index().drop(columns = 'index')\n        circle5 = df.query('xr < 0 & trakus_index > @trackus_threshold').reset_index().drop(columns = 'index')\n        chute['pp'] = chute['xr']\n        stretch2['pp'] = stretch2['xr']\n        circle3['pp'] = stretch+r*np.arctan2( circle3['xr'] - stretch, - circle3['yr'] ) \n        stretch4['pp'] = stretch + r*math.pi + (stretch - stretch4['xr'])       \n        circle5['pp'] = 2*stretch + r*math.pi + r*np.arctan2( - circle5['xr'] , circle5['yr'] ) \n\n        df = pd.concat([chute, stretch2, circle3, stretch4, circle5])\n        df = df.reset_index().drop(columns = ['index'])\n        return df , True\n    \n    else:\n        return df, False","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:20.177823Z","iopub.execute_input":"2022-11-09T11:44:20.178984Z","iopub.status.idle":"2022-11-09T11:44:20.191599Z","shell.execute_reply.started":"2022-11-09T11:44:20.178936Z","shell.execute_reply":"2022-11-09T11:44:20.190147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sprints = []\nfor i in tqdm(dict_t.keys()): \n    frame_now, result = pp_sprint(dict_t[i][0].query('distance_id < 705'), dict_t[i][1], dict_t[i][2])\n    sprints.append(frame_now)\nsprints_pp = pd.concat(sprints).reset_index(drop = True)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:20.203080Z","iopub.execute_input":"2022-11-09T11:44:20.203466Z","iopub.status.idle":"2022-11-09T11:44:27.535902Z","shell.execute_reply.started":"2022-11-09T11:44:20.203433Z","shell.execute_reply":"2022-11-09T11:44:27.534776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sprints_pp = sprints_pp.drop(columns=['latitude','longitude','x','y'])","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:27.537176Z","iopub.execute_input":"2022-11-09T11:44:27.537498Z","iopub.status.idle":"2022-11-09T11:44:27.945002Z","shell.execute_reply.started":"2022-11-09T11:44:27.537470Z","shell.execute_reply":"2022-11-09T11:44:27.943882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here I just check that all the chutes are having a shape that this model understands. \nThis plot is simply to make sure that indeed all the races I am looking at have horizontal chutes, which is confirmed. I hid it because it is here just for checking. ","metadata":{}},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (8,8)\nsprints_pp.plot(x='xr', y='yr', kind='scatter' )\nax = plt.gca()\nplt.gca().set_aspect('equal', adjustable='box')\nplt.show()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:27.946426Z","iopub.execute_input":"2022-11-09T11:44:27.947471Z","iopub.status.idle":"2022-11-09T11:44:33.234941Z","shell.execute_reply.started":"2022-11-09T11:44:27.947419Z","shell.execute_reply":"2022-11-09T11:44:33.233782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Add the fihishing line positions in coordinates of pp (projected position). \nI make a dictionary, dict_f with the finish line positions in \"pp\" (projected position) coordinates. For AQU D, AQU O, BEL D, SAR D, the finish line position in \"pp\" coordinates can be just calculated from the track information taking into account that our coordinates system is set in such a way that the 0 coordinate is at the beginning of the first stretch. \n> The finish line position in \"pp\" coordinates is \"length of stretch\" + \"length of turn\" + \"head of stretch to finish line\". As an example, for Belmont main track, we have in feet: (1097 + 843) + 2020 + 1097 = 5057 feet, that is, 1541.37m. \n\nFor the other tracks in the three venues the finish line positions have to be calculated a bit more careful by using the finish line position of the D track and the difference in the translation between the D and the current track. \n\nBelow are the track info that I used for the finish positions that can be calculated right away. ","metadata":{}},{"cell_type":"markdown","source":"AQUEDUCT Main Track: 1 1/8 Miles Around\n> - 1155.5' Head of Stretch to Finish Line\n> - 330' Finish Line to First Turn\n> - 470' Turn radius\n> - 1484.5' Length of Turn\n\nAQUEDUCT Inner Track (outer turf)\n> - 1 Mile Around\n> - 1174.9' Head of Stretch to Finish Line\n> - 332.6' Finish Line to First Turn\n> - 350' Turn Radius\n> - 1132.5 +' Length of Turn","metadata":{}},{"cell_type":"markdown","source":"BELMONT\n> - Main Track: 1 1/2 Miles Around\n> - 1097' Head of Stretch to Finish Line\n> - 843' Finish Line to First Turn\n> - 643' Turn Radius\n> - 2020' Length of Turn","metadata":{}},{"cell_type":"markdown","source":"SARATOGA\nMain Track: 1 1/8 Miles Around\n> - 1144' Head of Stretch to Finish Line\n> - 341' Finish Line to First Turn\n> - ~470' Turn Radius\n> - 1485' Length of Turn\n\n","metadata":{}},{"cell_type":"markdown","source":"\n","metadata":{}},{"cell_type":"code","source":"# these are the finish line positions in function of the projected position for each track. \n# They are extracted from the available information on the tracks. \ndict_f ={}\ndict_f['aqu_D'] = 1257.45\ndict_f['aqu_O'] = 1162.78\ndict_f['bel_D'] = 1541.37\ndict_f['sar_D'] = 1253.9472","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:33.237479Z","iopub.execute_input":"2022-11-09T11:44:33.238329Z","iopub.status.idle":"2022-11-09T11:44:33.244360Z","shell.execute_reply.started":"2022-11-09T11:44:33.238278Z","shell.execute_reply":"2022-11-09T11:44:33.243118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# here I calculate the final pp for belmont Dirt from our fit as opposed of calculating it from the measurements given. \n# we got 1538 as opposed to 1541.37, a 0.2% difference of 3.37m. \n2*dict_t['bel_D'][1] + dict_t['bel_D'][2]*math.pi - 256.9464","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:33.245715Z","iopub.execute_input":"2022-11-09T11:44:33.246202Z","iopub.status.idle":"2022-11-09T11:44:33.263422Z","shell.execute_reply.started":"2022-11-09T11:44:33.246166Z","shell.execute_reply":"2022-11-09T11:44:33.261990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_f['bel_T'] = 2*dict_t['bel_T'][1] + dict_t['bel_T'][2]*math.pi - 256.9464 - dict_t['bel_D'][3] + dict_t['bel_T'][3]\ndict_f['bel_I'] = 2*dict_t['bel_I'][1] + dict_t['bel_I'][2]*math.pi - 256.9464 - dict_t['bel_D'][3] + dict_t['bel_I'][3] \ndict_f['sar_T'] = 2*dict_t['sar_T'][1] + dict_t['sar_T'][2]*math.pi - 103.9368 - dict_t['sar_D'][3] + dict_t['sar_T'][3] \ndict_f['sar_I'] = 2*dict_t['sar_I'][1] + dict_t['sar_I'][2]*math.pi - 103.9368 - dict_t['sar_D'][3] + dict_t['sar_I'][3] ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:33.265013Z","iopub.execute_input":"2022-11-09T11:44:33.265453Z","iopub.status.idle":"2022-11-09T11:44:33.279557Z","shell.execute_reply.started":"2022-11-09T11:44:33.265420Z","shell.execute_reply":"2022-11-09T11:44:33.278258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# this is the dictionary with the pp values for the finishing line for all tracks and courses\ndict_f ","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:33.281554Z","iopub.execute_input":"2022-11-09T11:44:33.282058Z","iopub.status.idle":"2022-11-09T11:44:33.294250Z","shell.execute_reply.started":"2022-11-09T11:44:33.282012Z","shell.execute_reply":"2022-11-09T11:44:33.293120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Define a function to add actual speeds (sa) and projected speeds (sp) to each individual horse race\nI defined a function that adds the actual speed, 'sa' and the projected speed 'sp' to the dataframe corresponding to a horse racing. It smoothes the speed obtained from the data by applying a savgol filter (Mercier and Aftalion 2020). The speed can be added only to a dataframe containing the race of one horse. ","metadata":{}},{"cell_type":"code","source":"def add_vs(df):\n    '''\n    This function adds the actual speed, 'sa' and the projected speed 'sp' to the dataframe. It smoothes the speed obtained \n    from the data by applying a savgol filter as in the PLOS ONE paper. \n    Speeds can only be added to a single horse race\n    We add the sa, the actual speed, and sp the projected speed. the actual speed has to be larger than the projected speed\n    sp = sa cos(alpha) where alpha is the angle the instantaneous velocity makes with the direction paralel to the optimal \n    track direction\n    the vxs and xys , the actual speeds smoothed will be needed later when we calculate the draft.\n    '''\n    df = df.sort_values(by=['trakus_index']).reset_index().drop(columns = ['index'])\n    df_v = df[['xr','yr', 'pp','trakus_index']]\n    df_v = 4*df_v.diff()\n    df_v = df_v.rename(columns = {'xr': 'vx', 'yr': 'vy', 'pp': 'sp'})\n    df_v = df_v.drop(columns = ['trakus_index'])\n    df_v = df_v.fillna(0)\n    df = pd.merge(df, df_v, left_index=True, right_index=True)\n    df['sa'] = np.sqrt(df_v['vx']**2 + df_v['vy']**2)\n    sa = savgol_filter(df['sa'], 15, 3)\n    sp = savgol_filter(df['sp'], 15, 3)   \n    vxs = savgol_filter(df['vx'], 15, 3)\n    vys = savgol_filter(df['vy'], 15, 3)\n    \n    df['sas'] = sa\n    df['sps'] = sp\n    df['vxs'] = vxs\n    df['vys'] = vys\n    \n    if abs(df['sa'].max()) > 50 and abs(df['sa'].min()) > 50 :\n        df['sa'] = (df['sas']/df['sa'])*1000\n    \n    df['time'] = df['trakus_index']/4\n    if df['sa'].max() < 40:\n        status = True\n    else:\n        status = False\n    return df , status","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:33.296212Z","iopub.execute_input":"2022-11-09T11:44:33.296695Z","iopub.status.idle":"2022-11-09T11:44:33.309560Z","shell.execute_reply.started":"2022-11-09T11:44:33.296647Z","shell.execute_reply":"2022-11-09T11:44:33.308008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The plot below shows, for illustration, the actual speed and the projected speed (lower) for a horse chosen at random. \n\nNotice how a the turn the actual speed is higher than the \nprojected speed as expected. We have now the actual position, the projected position, and we can compare the total distance ran by a horse with the projected distance. We can work with actual velocity and with projected velocity.","metadata":{}},{"cell_type":"code","source":"#this is a rough way to select a one horse race\ncounter = 0\nfor element, dff in sprints_pp.groupby(['race_date', 'race_number', 'program_number']):\n    counter += 23\n    if counter == 14: \n        break\n\ndatenow = element[0]\nracenow = element[1]\nhorsenow = element[2]\ndff = sprints_pp.query('race_date == @datenow & race_number == @racenow & program_number == @horsenow')\ndff, state = add_vs(dff)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:33.311729Z","iopub.execute_input":"2022-11-09T11:44:33.312167Z","iopub.status.idle":"2022-11-09T11:44:34.873552Z","shell.execute_reply.started":"2022-11-09T11:44:33.312133Z","shell.execute_reply":"2022-11-09T11:44:34.872121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (12,12)\nplt.scatter(dff['time'], dff['sas'], color = 'r', label='actual speed')\nplt.scatter(dff['time'], dff['sps'], color = 'g', label='speed projected on the racetrack')\nplt.xlabel('time')\nplt.ylabel('speed (m/s)')\nplt.title('actual speed in red and projected speed in green for a random one horse trajectory')\nplt.legend(loc='center')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:34.875391Z","iopub.execute_input":"2022-11-09T11:44:34.876094Z","iopub.status.idle":"2022-11-09T11:44:35.107632Z","shell.execute_reply.started":"2022-11-09T11:44:34.876055Z","shell.execute_reply":"2022-11-09T11:44:35.106430Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Features definition and correlations with average velocity and position at finish","metadata":{}},{"cell_type":"markdown","source":"## We define a few quantities to be extracted from the data\n#### Velocity profile features\n> -  **delta_v_01** is difference between the maximum projected speed during the first stretch and the initial speed (usually 0)\n> -  **delta_v_21** is the difference between the projected speed at the end of the first turn and the projected speed at the end of the firt stretch\n> -  **delta\\_v_32** is the difference between the projected speed at the end of the race and the projected speed at the end of the turn\n> -  **delta_v_st** is the difference between maximum projected speed on the first stretch and the minimum projected speed on the first turn. It is an indication of the \"speed drop\" from the stretch to the turn\n> -  **delta_v_ts** is the difference between the maximum projected speed on the turn and the minimum projected speed on the last stretch. It is an indication of how well the horse can keep the distance with which it enters the last stretch\n> - **rel_sigma_v** is the standard deviation of the projected velocity divided by the average velocity for the race. It is an indication of how much the speed \"fluctuates\" from the average and according to our theoretical considerations the lower the **rel_sigma_v** the lower the power expenditure during the race for the same average speed. We expect low **rel_sigma_v** to correlate with lower numerical (better) position at finish and higher average speeds. \n\n#### Energy minimization features\n\n> - **draft percentage** is the percentage of the aerodynamic drag that is minimized during the race. We followed (Spence et al. 2012) in defining the draft position in the sense that we considered the horse as being \"shielded\" when it was at a distance less than 5m behind another horse and in a cone of $10^{o}$. We have not followed (Spence et al. 2012) in defining the draft percentage because their definition does not take into account the speed variation during the race.  The draft percentage I used takes into account that being shielded at low speeds is less impactful than being shielded at high speeds because the aerodynamic drag goes like $\\propto v^{3}$. The draft I am going to use is $draft_{v3} = \\frac {\\sum _{d} v_{d}^{3}}{\\sum _{i} v_{i}^{3}}$, where $i$ indicates the indices of all velocities and $d$ indicates only the indices of the velocities when the horse is in a so called \"drafting position\", that is less than 5m behind at least another horse at an angle no larger than $10^{o}$. The formula used here and the formula used in (Spence et al. 2012) give the same result at 0 and 1. \n> The draft percentage used in (Spence et al. 2012) is $draft_{old}= \\frac {\\sum _{d} 1}{\\sum _{i} 1}$ with the same definition for indices $i$ and $d$ as above.\n> - **delta_d** is the extra distance ran by the horse divided by the total projected distance of the race. It is an indication of how close to the rails the horse has been running during the turns to minimize total distance ran. Smaller is better. \n\n","metadata":{}},{"cell_type":"markdown","source":"# Which kind of features are more important in the outcome of an individual race? The outcomes are projected average speed (pa) and positon at finish.","metadata":{}},{"cell_type":"markdown","source":"I defined a function to add these features to each individual horse race. I will also add a \"draft\" feature, equal to the draft percentage as defined above. The result will be a dataframe containing all individual horse races and the extra features defined for each of these races. ","metadata":{}},{"cell_type":"code","source":"def add_info(df): # xc2 is the stretch length\n    row_list = []\n    listbad = []\n    for element, frame in tqdm(df.groupby(['track_id','race_date', 'race_number', 'program_number', \\\n        'course_type', 'track_condition', 'run_up_distance', 'race_type', 'purse', 'post_time', 'weight_carried', \\\n        'jockey', 'odds', 'position_at_finish', 'distance', 'win_prob', \\\n                                          \n                                          ])):\n        key = element[0].lower() + '_' + frame.course_type.unique()[0]\n        fp = dict_f[key]\n        xc2 = dict_t[key][1]\n        l = dict_t[key][2]*math.pi\n        \n        try:\n            frame, state = add_vs(frame)\n\n            # here i check wether the horse reached up to the finish line\n            end = frame.loc[frame.pp.transform(lambda z: (z >= fp) & (z.shift(1) < fp))]\n            if len(end) > 0:\n                df_race = frame.query('pp < @fp').sort_values(by=['trakus_index']).reset_index()    \n                sp = df_race['sps'].mean()\n                sa = df_race['sas'].mean()\n                delta_d = (df_race['sas'].sum()-df_race['sps'].sum())*100/df_race['sps'].sum()\n                distance_calc = df_race['pp'].max() - df_race['pp'].min()\n                sigma = df_race['sas'].std()\n                rel_sigma = sigma/sa\n\n                race_1st_stretch = df_race.query('pp < @xc2').sort_values(by=['trakus_index']).reset_index(drop=True)    \n                race_turn = df_race.query('pp > @xc2 & pp < @xc2 + @l').sort_values(by=['trakus_index']).reset_index(drop=True)    \n                race_2nd_stretch = df_race.query ('pp > @xc2 + @l & pp < @fp').sort_values(by=['trakus_index']).reset_index(drop=True)    \n                delta_v_st = race_1st_stretch['sps'].max() - race_turn['sps'].min()\n                delta_v_ts = race_turn['sps'].max() - race_2nd_stretch['sps'].min()\n                p1 = race_1st_stretch.tail(1).reset_index(drop=True)\n                p2 = race_turn.tail(1).reset_index(drop=True)\n                p3 = race_2nd_stretch.tail(1).reset_index(drop=True)\n\n                try:\n                    delta_v10 = p1.at[0,'sas'] -  race_1st_stretch['sps'].max()\n                except Exception:\n                    delta_v10 = None\n                try:    \n                    delta_v21 = p2.at[0,'sas'] - p1.at[0,'sas']\n                except Exception:\n                    delta_v21 = None\n                try:      \n                    delta_v32 = p3.at[0,'sas'] - p2.at[0, 'sas']\n                except Exception:\n                    delta_v32 = None\n                row = element + (sa,)\n                row = row + (sp,)\n                row = row + (delta_d,)\n                row = row + (distance_calc,)\n                row = row + (sigma,)\n                row = row + (rel_sigma,)\n                row = row + (delta_v_st,)\n                row = row + (delta_v_ts,)\n                row = row + (delta_v10,)\n                row = row + (delta_v21,)\n                row = row + (delta_v32,)\n                #row = row + (df_race.at[0,'position_at_finish'],)\n                row_list.append(row)\n            elif len(end) == 0:\n                listbad.append(element)\n\n        except Exception:\n            listbad.append(element)       \n        \n        \n        \n    df = pd.DataFrame(row_list, columns = ['track_id','race_date', 'race_number', 'program_number', \\\n        'course_type', 'track_condition', 'run_up_distance', 'race_type', 'purse', 'post_time', 'weight_carried', \\\n        'jockey', 'odds', 'position_at_finish', 'distance', 'win_prob', 'aa', 'pa', 'delta_d','distance_calc', 'sigma_v', \\\n        'rel_sigma_v', 'delta_v_st', 'delta_v_ts', 'delta_v10' , 'delta_v21', 'delta_v32'])\n                                              \n    return df  ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:35.109378Z","iopub.execute_input":"2022-11-09T11:44:35.110210Z","iopub.status.idle":"2022-11-09T11:44:35.133604Z","shell.execute_reply.started":"2022-11-09T11:44:35.110163Z","shell.execute_reply":"2022-11-09T11:44:35.132513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def draft_race(a_race):\n    '''\n    this function adds the v3s drafted to all the horses in a single race and returns the enriched dataframe\n    it does not return the draft as a number for that to happen another function has to be applied on this df\n    THIS FUNCTION IS FOR ONE RACE\n    It also adds the \"simple draft\"\n    '''\n    df_list = []\n    draft_list = []\n    speed_list = []\n    for element, df in a_race.groupby('program_number'):\n        \n        other_horses = set(list(a_race.program_number.unique())) - {element}\n        try:\n            df, state = add_vs(df)\n            merged = df.copy()\n            m_list = []\n            for o in other_horses:\n                df_o = a_race.query('program_number == @o')\n                df_o = df_o[['trakus_index','xr', 'yr']]\n                df_o =  df_o.rename(columns = {'xr': 'xr1', 'yr': 'yr1'})\n                merged = pd.merge(merged, df_o, on = \"trakus_index\")\n                merged['d_sqr'] = (merged['xr'] - merged['xr1'])**2 + (merged['yr'] - merged['yr1'])**2\n                merged['bool_d1'] = merged['d_sqr'] < 25\n                merged['dot'] = ((merged['xr1']-merged['xr'])*merged['vxs']+(merged['yr1']-merged['yr'])*merged['vys']) / (np.sqrt((merged['xr1']-merged['xr'])**2 + (merged['yr1']-merged['yr'])**2)*np.sqrt(merged['vxs']**2 + merged['vys']**2))\n                limit = math.cos(3.1415*10/180)\n                merged['bool_angle'] = merged['dot'] > limit\n                merged['draft' + str(o)] = (merged['bool_d1']*merged['bool_angle'])\n                merged = merged.drop( columns = ['d_sqr', 'dot', 'bool_angle', 'xr1', 'yr1'])\n                m_list.append('draft' + str(o))\n            merged['draft'] = merged[m_list].sum(axis = 1)\n            merged['draft'] = merged['draft'].astype(bool)\n            #merged['draft'] = merged['draft']*merged['bool_d1']\n            #merged = merged.drop( columns = ['bool_d1'])\n            #merged = merged.drop(columns = m_list)\n            \n            merged['v3'] = merged['sas']**3\n            merged['draft_v3'] = merged['v3']*merged['draft']\n            merged['draft_simple'] = merged['sas']*merged['draft']     \n            the_draft = merged['draft_v3'].sum()/merged['v3'].sum()\n            \n            df_list.append(merged)\n        except Exception:\n            pass\n    #print(m_list)\n    \n    final_df = pd.concat(df_list)\n\n    return final_df    \n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:35.134955Z","iopub.execute_input":"2022-11-09T11:44:35.136043Z","iopub.status.idle":"2022-11-09T11:44:35.153654Z","shell.execute_reply.started":"2022-11-09T11:44:35.136009Z","shell.execute_reply":"2022-11-09T11:44:35.152491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def draft_all(df1):\n    '''\n    this function adds the velocities at the third power v^3 only when the horse is in draft position. to find the numerical\n    draft an extra operation of summing is needed.\n    it adds them to the whole dataframe , more horses such as aqu_sprint_d\n    '''\n    dfs = []\n    for element, df in tqdm(df1.groupby(['track_id','race_date', 'race_number'])):\n        df_race = draft_race(df)\n        dfs.append(df_race)\n    return pd.concat(dfs)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:35.154961Z","iopub.execute_input":"2022-11-09T11:44:35.155956Z","iopub.status.idle":"2022-11-09T11:44:35.169457Z","shell.execute_reply.started":"2022-11-09T11:44:35.155922Z","shell.execute_reply":"2022-11-09T11:44:35.168053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_draft(df_all):\n    '''\n    this function simply sums up the v3 draft and v3s from the dataframe that came out of add_draftv3 and returns a dataframe \n    with the race time, number program number and numerical draft. \n    '''\n    row_list = []\n    list_bad = []\n    \n    for element, frame in tqdm(df_all.groupby(['track_id','race_date', 'race_number', 'program_number','course_type'])):\n        key = element[0].lower()+'_'+element[4]\n        fp = dict_f[key]\n        \n        end = frame.loc[frame.pp.transform(lambda z: (z >= fp) & (z.shift(1) < fp))]\n        if end is not None:\n            df_race = frame.query('pp < @fp').sort_values(by=['trakus_index']).reset_index()     \n            the_draft = df_race['draft'].sum()/len(df_race)\n            draft_v3  = df_race['draft_v3'].sum()/df_race['v3'].sum()\n            row = element + (the_draft,)    \n            row = row +(draft_v3, )\n            row_list.append(row)\n        elif end is None:\n            listbad.append(element)\n\n        \n    df = pd.DataFrame(row_list, columns = ['track_id','race_date', 'race_number', 'program_number', 'course_type','draft_simple','draft_v3'])\n    \n    return df \n    ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:35.170878Z","iopub.execute_input":"2022-11-09T11:44:35.171256Z","iopub.status.idle":"2022-11-09T11:44:35.181670Z","shell.execute_reply.started":"2022-11-09T11:44:35.171223Z","shell.execute_reply":"2022-11-09T11:44:35.180817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# here I add all the extra information to each individual horse race\ndf_corr = add_info(sprints_pp)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:44:35.183051Z","iopub.execute_input":"2022-11-09T11:44:35.183816Z","iopub.status.idle":"2022-11-09T11:49:45.583051Z","shell.execute_reply.started":"2022-11-09T11:44:35.183771Z","shell.execute_reply":"2022-11-09T11:49:45.581917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_d(df):\n    listdf = []\n    for element, frame in tqdm(df.groupby(['track_id','race_date', 'race_number'])):\n        frame['delta_d_race'] = (frame['delta_d'] - frame['delta_d'].mean())/frame['delta_d'].mean()\n        listdf.append(frame)\n    df = pd.concat(listdf)\n    df = df.reset_index(drop=True)\n    return df","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:49:45.584552Z","iopub.execute_input":"2022-11-09T11:49:45.584951Z","iopub.status.idle":"2022-11-09T11:49:45.591111Z","shell.execute_reply.started":"2022-11-09T11:49:45.584918Z","shell.execute_reply":"2022-11-09T11:49:45.589874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr = add_d(df_corr)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:49:45.593001Z","iopub.execute_input":"2022-11-09T11:49:45.593488Z","iopub.status.idle":"2022-11-09T11:49:47.309238Z","shell.execute_reply.started":"2022-11-09T11:49:45.593444Z","shell.execute_reply":"2022-11-09T11:49:47.308095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr = df_corr.dropna()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:49:47.311109Z","iopub.execute_input":"2022-11-09T11:49:47.311575Z","iopub.status.idle":"2022-11-09T11:49:47.324714Z","shell.execute_reply.started":"2022-11-09T11:49:47.311527Z","shell.execute_reply":"2022-11-09T11:49:47.323544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sprints_draft1 = draft_all(sprints_pp)\nsprints_draft2 =add_draft(sprints_draft1)\ndf_draft = sprints_draft2.copy()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T11:49:47.325979Z","iopub.execute_input":"2022-11-09T11:49:47.327074Z","iopub.status.idle":"2022-11-09T12:08:32.171192Z","shell.execute_reply.started":"2022-11-09T11:49:47.327037Z","shell.execute_reply":"2022-11-09T12:08:32.170261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_corr = pd.merge(df_corr, df_draft,  how= 'left', left_on= ['race_date', 'race_number', 'program_number','track_id','course_type'], \\\nright_on = ['race_date', 'race_number', 'program_number','track_id','course_type'] )\n\ndf_corr2 = final_corr[[ 'weight_carried', 'position_at_finish', 'distance',\n       'pa', 'delta_d', 'sigma_v', 'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts', 'delta_v10', 'delta_v21',\n       'delta_v32', 'delta_d_race', 'draft_simple', 'draft_v3'\n    ]].copy()\n\n\ndf_corr2  = df_corr2[(np.abs(stats.zscore(df_corr2)) < 3).all(axis=1)] ","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:32.172832Z","iopub.execute_input":"2022-11-09T12:08:32.173223Z","iopub.status.idle":"2022-11-09T12:08:32.214267Z","shell.execute_reply.started":"2022-11-09T12:08:32.173187Z","shell.execute_reply":"2022-11-09T12:08:32.213030Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The dataframe containing all the individual horse races info is called final_corr. Below are the data types for each column. ","metadata":{}},{"cell_type":"code","source":"final_corr.dtypes","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:32.216035Z","iopub.execute_input":"2022-11-09T12:08:32.216403Z","iopub.status.idle":"2022-11-09T12:08:32.225413Z","shell.execute_reply.started":"2022-11-09T12:08:32.216371Z","shell.execute_reply":"2022-11-09T12:08:32.224326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We apply groupby to see how many races (that is, an actual race ran by a few horses) we have. We have 1134 out of 2000 races, which is good, taking into account that we are looking only at short races, out of three possible types of races: short, medium, and long.","metadata":{}},{"cell_type":"code","source":"df_corr.groupby(['track_id','race_date','race_number']).count()","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:32.226919Z","iopub.execute_input":"2022-11-09T12:08:32.227278Z","iopub.status.idle":"2022-11-09T12:08:32.288341Z","shell.execute_reply.started":"2022-11-09T12:08:32.227241Z","shell.execute_reply":"2022-11-09T12:08:32.287031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We start with correlation plots which are a first look on where the effects on position at finish and average speed may lie. Because there are many features we split them in two plots, keeping **pa** and **position_at_finish** in both of them.\n\nHere's the first correlations plot with features that do not correlate well with either the average projected speed (**pa**) or the **position_at_finish**.","metadata":{}},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (8,8)\n\nsns.pairplot(data=df_corr2, vars = [ 'weight_carried',  'distance',\n       'delta_d', 'delta_v_ts', 'delta_v10', \n       'delta_v32',  'pa', 'position_at_finish'  ] )\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:32.290713Z","iopub.execute_input":"2022-11-09T12:08:32.291200Z","iopub.status.idle":"2022-11-09T12:08:45.558157Z","shell.execute_reply.started":"2022-11-09T12:08:32.291150Z","shell.execute_reply":"2022-11-09T12:08:45.556990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here's the heatmap corresponding to the correlation plot above, where we can read the correlation coefficients. ","metadata":{}},{"cell_type":"code","source":"df_corrs = df_corr2[[  'weight_carried',  'distance',\n       'delta_d', 'delta_v_ts', 'delta_v10', \n       'delta_v32',  'pa', 'position_at_finish' ]].corr()\nsns.heatmap(df_corrs, annot=True)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:45.559675Z","iopub.execute_input":"2022-11-09T12:08:45.560013Z","iopub.status.idle":"2022-11-09T12:08:46.146913Z","shell.execute_reply.started":"2022-11-09T12:08:45.559982Z","shell.execute_reply":"2022-11-09T12:08:46.145634Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The second correlations plot contains features that correlate well both with the average projected speed (**pa**) and the **position_at_finish**. ","metadata":{}},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (8,8)\n\nsns.pairplot(data=df_corr2, vars = [ 'rel_sigma_v',\n       'delta_v_st', 'delta_v21',\n        'draft_v3','position_at_finish', \n       'pa'] )\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:46.148505Z","iopub.execute_input":"2022-11-09T12:08:46.148884Z","iopub.status.idle":"2022-11-09T12:08:53.940945Z","shell.execute_reply.started":"2022-11-09T12:08:46.148852Z","shell.execute_reply":"2022-11-09T12:08:53.939391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corrs = df_corr2[[   'rel_sigma_v',\n       'delta_v_st', 'delta_v21',\n        'draft_v3','position_at_finish', \n       'pa' ]].corr()\nsns.heatmap(df_corrs, annot=True)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:53.942484Z","iopub.execute_input":"2022-11-09T12:08:53.942896Z","iopub.status.idle":"2022-11-09T12:08:54.376340Z","shell.execute_reply.started":"2022-11-09T12:08:53.942860Z","shell.execute_reply":"2022-11-09T12:08:54.375169Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### First observations from the heatmaps above:\n> - The position at finish correlates well with the average projected speed, **pa**, the higher the ** pa ** the lower the position at finish, which is to be expected. \n> - The position at finish does not correlate well with the draft percentage\n> - the **delta_d** doesn't seem to help much either, while it correlates negatively with the average projected speed **pa** and positively with the position at finish (where small is good) the correlations are really small. \n> - The speed drops, **delta_v_st** and **delta_v_ts** strongly inversely correlate  with **pa** and correlate with **position_at_finish** (larger position at finish). \n> - **delta_v10**, **delta_v21**, and **delta_v32** also correlate well with **pa** and inversely correlate with **positon_st_finish** (smaller position at finish).\n\nThe general take home message is that \"strategic racing\", that is sitting long in a drafting position or keeping the total distance ran as short as possible does not seem to work and **it's not worth the risk of injury**.","metadata":{}},{"cell_type":"markdown","source":"Indeed the higher the average projected speed the lower position at finish so we can also use the average speed as an indication \nof a racing strategy success. Notice that the error bars are quite large, as opposed to the error bars in (Spence et al. 2012).","metadata":{}},{"cell_type":"code","source":"list_p = []\nlist_ap = []\nlist_ap_std = []\nfor i, df in df_corr2.groupby('position_at_finish'):\n    average = df['pa'].mean()\n    std_dev = df['pa'].std()\n    list_p.append(i)\n    list_ap.append(average)\n    list_ap_std.append(std_dev)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.377791Z","iopub.execute_input":"2022-11-09T12:08:54.378125Z","iopub.status.idle":"2022-11-09T12:08:54.392902Z","shell.execute_reply.started":"2022-11-09T12:08:54.378094Z","shell.execute_reply":"2022-11-09T12:08:54.391694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (8,8)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.394704Z","iopub.execute_input":"2022-11-09T12:08:54.395049Z","iopub.status.idle":"2022-11-09T12:08:54.399913Z","shell.execute_reply.started":"2022-11-09T12:08:54.395018Z","shell.execute_reply":"2022-11-09T12:08:54.398928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.scatter(x=list_p, y=list_ap ,marker = \"o\",  linewidths = 1, color = 'g')\nplt.errorbar(list_p, list_ap, yerr=list_ap_std, fmt=\"o\")\nplt.xlabel('position at finish')\nplt.ylabel('average speed')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.401294Z","iopub.execute_input":"2022-11-09T12:08:54.401601Z","iopub.status.idle":"2022-11-09T12:08:54.625959Z","shell.execute_reply.started":"2022-11-09T12:08:54.401573Z","shell.execute_reply":"2022-11-09T12:08:54.624665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our data on drafting does not corroborate previously reported data (Spence et al. 2012). Drafting seems to have no effect on the position at finish. ","metadata":{}},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (8,8)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.627472Z","iopub.execute_input":"2022-11-09T12:08:54.628927Z","iopub.status.idle":"2022-11-09T12:08:54.634486Z","shell.execute_reply.started":"2022-11-09T12:08:54.628878Z","shell.execute_reply":"2022-11-09T12:08:54.633293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_p = []\nlist_ap = []\nlist_ap_std = []\nfor i, df in df_corr2.groupby('position_at_finish'):\n    average = df['draft_v3'].mean()\n    std_dev = df['draft_v3'].std()\n    list_p.append(i)\n    list_ap.append(average)\n    list_ap_std.append(std_dev)\n\nplt.scatter(x=list_p, y=list_ap ,marker = \"o\",  linewidths = 1, color = 'g')\nplt.errorbar(list_p, list_ap, yerr=list_ap_std, fmt=\"o\")\nplt.xlabel('position at finish')\nplt.ylabel('the draft')\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.636659Z","iopub.execute_input":"2022-11-09T12:08:54.637382Z","iopub.status.idle":"2022-11-09T12:08:54.873357Z","shell.execute_reply.started":"2022-11-09T12:08:54.637343Z","shell.execute_reply":"2022-11-09T12:08:54.872027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Gradient boosted trees confirm the preliminary findings and add more insight","metadata":{}},{"cell_type":"markdown","source":"Gradient boosted trees are a great tool to identify the features that exert the largest influence on the outcome, either average projected speed or position at finish. Although, as you will see they can be used to predict the race outcome, by using the SHAP values we will determine which features exert the most influence on the race outcome and in which way. This information is valuable for us to determine which type of racing strategies are to be pursued. The advantage of XGBoost is the fact that it can also take into account the effect of categorical features such as course_type and track_condition. \nWe are going to first do a regression fit on the **ordinal** feature that is the position at finish. ","metadata":{}},{"cell_type":"markdown","source":"#### Here I eliminate the outliers from the numerical values\nThe jockey names will not enter the models because it will mean creating a very large number of categorical data features. The outliers from the numerical values are eliminated. The resulting dataframe is called final_corr. ","metadata":{}},{"cell_type":"code","source":"final_corr_s = final_corr[['track_id','race_date','program_number','course_type','track_condition', 'race_type']]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.874997Z","iopub.execute_input":"2022-11-09T12:08:54.875364Z","iopub.status.idle":"2022-11-09T12:08:54.883099Z","shell.execute_reply.started":"2022-11-09T12:08:54.875330Z","shell.execute_reply":"2022-11-09T12:08:54.881533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_corr_n = final_corr[[ 'race_number', \n       'run_up_distance', \n       'purse', 'post_time', 'weight_carried',  'odds',\n       'position_at_finish', 'distance', 'win_prob', 'aa', 'pa', 'delta_d',\n       'distance_calc', 'sigma_v', 'rel_sigma_v', 'delta_v_st', 'delta_v_ts',\n       'delta_v10', 'delta_v21', 'delta_v32', 'delta_d_race', 'draft_simple',\n       'draft_v3']]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.884594Z","iopub.execute_input":"2022-11-09T12:08:54.885048Z","iopub.status.idle":"2022-11-09T12:08:54.898231Z","shell.execute_reply.started":"2022-11-09T12:08:54.885006Z","shell.execute_reply":"2022-11-09T12:08:54.897117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"set(final_corr.columns) - set(final_corr_s.columns) - set(final_corr_n.columns)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.899731Z","iopub.execute_input":"2022-11-09T12:08:54.900687Z","iopub.status.idle":"2022-11-09T12:08:54.914541Z","shell.execute_reply.started":"2022-11-09T12:08:54.900649Z","shell.execute_reply":"2022-11-09T12:08:54.913305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_corr_n = final_corr_n[(np.abs(stats.zscore(final_corr_n)) < 3).all(axis=1)]\nfinal_corr = pd.merge(final_corr_n,final_corr_s,left_index=True, right_index=True)\nfinal_corr['furlongs'] = final_corr['distance']/100","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:08:54.915950Z","iopub.execute_input":"2022-11-09T12:08:54.916310Z","iopub.status.idle":"2022-11-09T12:08:54.939789Z","shell.execute_reply.started":"2022-11-09T12:08:54.916273Z","shell.execute_reply":"2022-11-09T12:08:54.938690Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Let's model the average projected speed **pa** on both categorical and numerical features ","metadata":{"execution":{"iopub.status.busy":"2022-11-03T13:31:41.177232Z","iopub.execute_input":"2022-11-03T13:31:41.177774Z","iopub.status.idle":"2022-11-03T13:31:41.183023Z","shell.execute_reply.started":"2022-11-03T13:31:41.177732Z","shell.execute_reply":"2022-11-03T13:31:41.181954Z"}}},{"cell_type":"markdown","source":"The approach used here follows closely  https://www.kaggle.com/code/carlmcbrideellis/an-introduction-to-xgboost-regression","metadata":{}},{"cell_type":"code","source":"X = final_corr[['course_type', 'track_condition', 'race_type', 'purse','weight_carried', 'win_prob', 'distance', \\\n'delta_d', 'rel_sigma_v', 'delta_v_st', 'delta_v_ts', 'delta_v10', 'delta_v21', 'delta_v32', 'draft_v3' \\\n                 \n                 ]]\nX.index.name = 'index'","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:54.941455Z","iopub.execute_input":"2022-11-09T12:08:54.941807Z","iopub.status.idle":"2022-11-09T12:08:54.949889Z","shell.execute_reply.started":"2022-11-09T12:08:54.941775Z","shell.execute_reply":"2022-11-09T12:08:54.948510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = final_corr['pa']\ny.index.name = 'index'","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:54.964370Z","iopub.execute_input":"2022-11-09T12:08:54.964704Z","iopub.status.idle":"2022-11-09T12:08:54.977501Z","shell.execute_reply.started":"2022-11-09T12:08:54.964662Z","shell.execute_reply":"2022-11-09T12:08:54.976215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat_columns = ['course_type', 'track_condition', 'race_type']\nX = pd.get_dummies(X, columns = cat_columns )","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:54.978793Z","iopub.execute_input":"2022-11-09T12:08:54.979134Z","iopub.status.idle":"2022-11-09T12:08:54.998334Z","shell.execute_reply.started":"2022-11-09T12:08:54.979104Z","shell.execute_reply":"2022-11-09T12:08:54.996901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.2, random_state=7)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:55.000151Z","iopub.execute_input":"2022-11-09T12:08:55.000524Z","iopub.status.idle":"2022-11-09T12:08:55.011364Z","shell.execute_reply.started":"2022-11-09T12:08:55.000487Z","shell.execute_reply":"2022-11-09T12:08:55.010166Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"regressor=xgb.XGBRegressor(eval_metric='rmsle')","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:55.012930Z","iopub.execute_input":"2022-11-09T12:08:55.013333Z","iopub.status.idle":"2022-11-09T12:08:55.021559Z","shell.execute_reply.started":"2022-11-09T12:08:55.013296Z","shell.execute_reply":"2022-11-09T12:08:55.020427Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"param_grid = {\"max_depth\":    [4, 5, 6],\n              \"n_estimators\": [500, 600, 700],\n              \"learning_rate\": [0.01, 0.015, 0.02]}","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:55.023115Z","iopub.execute_input":"2022-11-09T12:08:55.023441Z","iopub.status.idle":"2022-11-09T12:08:55.034140Z","shell.execute_reply.started":"2022-11-09T12:08:55.023413Z","shell.execute_reply":"2022-11-09T12:08:55.032927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"search = GridSearchCV(regressor, param_grid, cv=5).fit(X_train, y_train)\n\nprint(\"The best hyperparameters are \",search.best_params_)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:08:55.035926Z","iopub.execute_input":"2022-11-09T12:08:55.036488Z","iopub.status.idle":"2022-11-09T12:20:51.200073Z","shell.execute_reply.started":"2022-11-09T12:08:55.036441Z","shell.execute_reply":"2022-11-09T12:20:51.198162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"regressor=xgb.XGBRegressor(learning_rate = search.best_params_[\"learning_rate\"],\n                           n_estimators  = search.best_params_[\"n_estimators\"],\n                           max_depth     = search.best_params_[\"max_depth\"],\n                           eval_metric='rmsle')\n\nregressor.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:20:51.202803Z","iopub.execute_input":"2022-11-09T12:20:51.204024Z","iopub.status.idle":"2022-11-09T12:20:59.896149Z","shell.execute_reply.started":"2022-11-09T12:20:51.203972Z","shell.execute_reply":"2022-11-09T12:20:59.894873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions = regressor.predict(X_test)\nmae_xgb = mean_absolute_error(y_test, predictions)\nprint(f'The model mean absolute error is {mae_xgb:0.2f} ')","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:20:59.897873Z","iopub.execute_input":"2022-11-09T12:20:59.898592Z","iopub.status.idle":"2022-11-09T12:20:59.923910Z","shell.execute_reply.started":"2022-11-09T12:20:59.898546Z","shell.execute_reply":"2022-11-09T12:20:59.922860Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I will look later how the XGBoost prediction compares with the baseline but until now, a mean absolute error of 0.21 is not bad for integer ordinal targets. For now, we just want to look now at the influence of various features on the fit. ","metadata":{}},{"cell_type":"markdown","source":"### Next step is to use the SHAP values for our features to see which features influence the **pa** the most ","metadata":{}},{"cell_type":"code","source":"explainer = shap.TreeExplainer(regressor)\nshap_values = explainer.shap_values(X_test)\nshap.summary_plot(shap_values, X_test, show = False)\nplt.title(\"SHAP values for the features influencing the average projectes speed pa\")\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:20:59.937124Z","iopub.execute_input":"2022-11-09T12:20:59.937467Z","iopub.status.idle":"2022-11-09T12:21:06.994568Z","shell.execute_reply.started":"2022-11-09T12:20:59.937436Z","shell.execute_reply":"2022-11-09T12:21:06.993451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### How to read the SHAP charts: \n> -  the size of the SHAP values , in line with the position in the SHAP plot, shows the **magnitude** of the influence that particular feature has on the model\n> - the color distribution of the SHAP value indicates whether the correlation is proportional or inversely proportional. Red means that a larger value for a particular feature means a larger value for the average projectes speed, blue is the other way round. To put it in simpler terms, red on the right means that feature moves in the same direction with the average projectes speed, blue on the right means that feature moves in opposite direction. ","metadata":{}},{"cell_type":"markdown","source":"#### Learnings from the SHAP values for the XGB model to predict **pa**, the Projected Average speed\n> - the larger the purse money, the higher the average speed. This observation totally makes sense because there will be more prize money for high performance horses\n> - a high **rel_sigma_v** means lower average projected speeds, in line with our theoretical considerations and the numerical findings in Annex II. In other words, the horse and jockey must try to minimize speed fluctuations as much as possible. It is nice when numerical and experimental data fit our theoretical considerations. \n> -  the **course types** and **track conditions** are of utmost importance in the average speed outcome, indicating that the approximation used in previous work (Spence et al. 2012) to put them in the model as random variables is not valid.\n> - sitting in a \"drafting position\" ** is not helping ** , finding that corroborates our previous findings from the correlation plots and **finish position** vs **draft**. \n> - **delta_10** correlates positively with the average speed, that is the highest speed upon the initial sprint is an important indication for the average race speed. \n> - the race distance expressed meters is a late 7th in importance regarding the average race speed, again in contradiction to the prevous findings (Spence et al. 2012), where the race distance was the most important model parameter. \n> - a number of features related to the speed profile during the race $s = s(d)$, such as **delta_v_st**, **delta_v_ts**, **delta_v_32** follow closely\n> - as shown in the correlations graph above, the path length, the ** delta_d** feature has very little impact on the average projected speed. ","metadata":{}},{"cell_type":"markdown","source":"### Let's model the **position at finish** on both categorical and numerical features \nModeling the final position requires selecting a subset of races and eliminating a number of categorical features. For example, it makes no sense to leave in the model the **course_type** and **track_condition** because these features are identical for all horses competing agains each other in a race. The only categorical feature left is the **race_type**. I will look at course_type D and track_condition FT simply because these races occur most often.","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"races_d_ft = final_corr.query('course_type == \"D\" & track_condition == \"FT \"')","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:21:06.999902Z","iopub.execute_input":"2022-11-09T12:21:07.000475Z","iopub.status.idle":"2022-11-09T12:21:07.019546Z","shell.execute_reply.started":"2022-11-09T12:21:07.000423Z","shell.execute_reply":"2022-11-09T12:21:07.018458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nr_horse_races = len(races_d_ft)\nnr_races = len(races_d_ft.groupby(['track_id','race_date','race_number']).count())\nprint(f'The total number of individual horse races is {nr_horse_races:g}, while the total number of races is {nr_races:g}' )","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:21:07.020954Z","iopub.execute_input":"2022-11-09T12:21:07.021804Z","iopub.status.idle":"2022-11-09T12:21:07.037929Z","shell.execute_reply.started":"2022-11-09T12:21:07.021762Z","shell.execute_reply":"2022-11-09T12:21:07.036977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X5 = races_d_ft[['race_type', 'purse','weight_carried', 'distance','win_prob', \\\n'delta_d', 'rel_sigma_v', 'delta_v_st', 'delta_v_ts', 'delta_v10', 'delta_v21', 'delta_v32', 'draft_v3' \\\n                 ]]\nX5.index.name = 'index'","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:21:07.039225Z","iopub.execute_input":"2022-11-09T12:21:07.040055Z","iopub.status.idle":"2022-11-09T12:21:07.048313Z","shell.execute_reply.started":"2022-11-09T12:21:07.040012Z","shell.execute_reply":"2022-11-09T12:21:07.046841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat_columns = [ 'race_type']\nX5 = pd.get_dummies(X5, columns = cat_columns )\ny5 = races_d_ft['position_at_finish']\ny5.index.name = 'index'","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:21:07.049911Z","iopub.execute_input":"2022-11-09T12:21:07.050322Z","iopub.status.idle":"2022-11-09T12:21:07.072018Z","shell.execute_reply.started":"2022-11-09T12:21:07.050285Z","shell.execute_reply":"2022-11-09T12:21:07.070488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train5, X_test5, y_train5, y_test5 = train_test_split(X5,y5, test_size = 0.2, random_state=7)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:21:07.073728Z","iopub.execute_input":"2022-11-09T12:21:07.074082Z","iopub.status.idle":"2022-11-09T12:21:07.084774Z","shell.execute_reply.started":"2022-11-09T12:21:07.074050Z","shell.execute_reply":"2022-11-09T12:21:07.083400Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"regressor=xgb.XGBRegressor(eval_metric='rmsle')\nparam_grid = {\"max_depth\":    [4, 5, 6],\n              \"n_estimators\": [500, 600, 700],\n              \"learning_rate\": [0.01, 0.015, 0.02]}\nsearch = GridSearchCV(regressor, param_grid, cv=5).fit(X_train5, y_train5)\n\nprint(\"The best hyperparameters are \",search.best_params_)\nregressor=xgb.XGBRegressor(learning_rate = search.best_params_[\"learning_rate\"],\n                           n_estimators  = search.best_params_[\"n_estimators\"],\n                           max_depth     = search.best_params_[\"max_depth\"],\n                           eval_metric='rmsle')\nregressor.fit(X_train5, y_train5)\n\npredictions = regressor.predict(X_test5)\n\nprint('mae is ' , mean_absolute_error(y_test5, predictions))\n\nprint( 'std ground truth is ', y_test5.std())","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:21:07.087073Z","iopub.execute_input":"2022-11-09T12:21:07.087741Z","iopub.status.idle":"2022-11-09T12:29:03.311218Z","shell.execute_reply.started":"2022-11-09T12:21:07.087652Z","shell.execute_reply":"2022-11-09T12:29:03.310173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Next step is to use the SHAP values for out features to see which features influence the **position_at_finish** the most ","metadata":{}},{"cell_type":"code","source":"explainer = shap.TreeExplainer(regressor)\nshap_values = explainer.shap_values(X_test5)\nshap.summary_plot(shap_values, X_test5, show = False)\nplt.title(\"SHAP values for the features influencing the ORDINAL feature, position_at_finish\")\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:29:03.312838Z","iopub.execute_input":"2022-11-09T12:29:03.313886Z","iopub.status.idle":"2022-11-09T12:29:06.671451Z","shell.execute_reply.started":"2022-11-09T12:29:03.313848Z","shell.execute_reply":"2022-11-09T12:29:06.670509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Because a large positon at finish means lower performance, this time features that \"help\" the final position have blue on the right side (positive values) and those that impede it have red on the right side of the graph.\n> - As shown before, the draft does not help performace, corroborating our previous findings\n> - The **win probability** correlates positively with a low numerical position at finish (good position), which is to be expected and only validates that our model is working properly\n> - The next features are the five defined velocity differences and they strongly correlate with the outcome. The **delta_v_st** and **delta_v_ts** correlation high numerical (bad) position at finish confirm that large velocity drops from stretch to turn and from turn to stretch predict low performance, and the forward looking **delta_v10**, **delta_v_21** and **delta_v_32** confirm that keeping a steady or increasing speed during the race predicts a good outcome. \n> - **delta_d** againd does not have a large predicting values, as we saw in our correlation plots and previous SHAP values plots. ","metadata":{}},{"cell_type":"markdown","source":"# We can beat the odds. That is, we can predict the position at finish better than the winning probabilities.\nThe most straightforward way to predict the position at finish is to simply rank the winning probability for each horse within a race and use that ranking as the position at finish value. This approach is the benchmark for any algorithm predicting the ** position at finish** outcome.<br>\nTo test wether our algorithm performs better than the naive \"odds algorithm\" described above, we run the predictor and rank each horse in the order of its predicted position at finish and then use those indices as the position at finish values. We do it again on course_type D and track_condition FT <br>\n","metadata":{}},{"cell_type":"markdown","source":"Now I make sure that in the training data and test data there are different races. To make sure that the races are separated between the training and testing data I make a list of all races and shuffle it and then select training races and testing races. After I select training races and testing races I put back together the dataframes for my XGB model. ","metadata":{}},{"cell_type":"code","source":"races_list = []\nfor element, df in races_d_ft.groupby(['track_id','race_date','race_number']):\n    races_list.append(element)\nrandom.shuffle(races_list)\nraces_list[2]","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:29:06.685905Z","iopub.execute_input":"2022-11-09T12:29:06.686741Z","iopub.status.idle":"2022-11-09T12:29:06.709562Z","shell.execute_reply.started":"2022-11-09T12:29:06.686695Z","shell.execute_reply":"2022-11-09T12:29:06.708460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_races = races_list[:int(0.8*len(races_list))]\ntest_races = races_list[int(0.8*len(races_list)):]","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:29:06.711465Z","iopub.execute_input":"2022-11-09T12:29:06.712218Z","iopub.status.idle":"2022-11-09T12:29:06.718564Z","shell.execute_reply.started":"2022-11-09T12:29:06.712179Z","shell.execute_reply":"2022-11-09T12:29:06.717214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_list = []\ntest_list = []\nfor element in train_races:\n    track_now = element[0]\n    date_now = element[1]\n    number_now = element[2]\n    train_list.append(races_d_ft.query('track_id == @track_now & race_date == @date_now & race_number == @number_now'))\nfor element in test_races:\n    track_now = element[0]\n    date_now = element[1]\n    number_now = element[2]\n    test_list.append(races_d_ft.query('track_id == @track_now & race_date == @date_now & race_number == @number_now'))    \ntrain_df = pd.concat(train_list)\ntest_df = pd.concat(test_list)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:29:06.720770Z","iopub.execute_input":"2022-11-09T12:29:06.721785Z","iopub.status.idle":"2022-11-09T12:29:09.554849Z","shell.execute_reply.started":"2022-11-09T12:29:06.721735Z","shell.execute_reply":"2022-11-09T12:29:09.553621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To train the model we first use ONLY features that describe variations BEFORE the race: 'weight_carried', 'win_prob', 'delta_v_st',  'delta_v10', 'delta_v21'. This actually means predicting the outcome at the end of the turn. Can we \"beat the odds\" in the literal sense, that is offer a better prediction than the odds?","metadata":{}},{"cell_type":"code","source":"X_train9 = train_df[['weight_carried', 'win_prob', 'delta_v_st',  'delta_v10', 'delta_v21']]\nX_train9.index.name = 'index'\n\ny_train9 = train_df['position_at_finish']\ny_train9.index.name = 'index'\n\ny_test9 = test_df[['position_at_finish']]\ny_test9.index.name ='index'\n\nX_test9 = test_df[['weight_carried', 'win_prob', 'delta_v_st',  'delta_v10', 'delta_v21']]\nX_test9.index.name = 'index'","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:29:09.578892Z","iopub.execute_input":"2022-11-09T12:29:09.579535Z","iopub.status.idle":"2022-11-09T12:29:09.592816Z","shell.execute_reply.started":"2022-11-09T12:29:09.579477Z","shell.execute_reply":"2022-11-09T12:29:09.591345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"regressor=xgb.XGBRegressor(eval_metric='rmsle')\nparam_grid = {\"max_depth\":    [4, 5, 6],\n              \"n_estimators\": [500, 600, 700],\n              \"learning_rate\": [0.01, 0.015, 0.02]}\nsearch = GridSearchCV(regressor, param_grid, cv=5).fit(X_train5, y_train5)\n\nprint(\"The best hyperparameters are \",search.best_params_)\n\nregressor=xgb.XGBRegressor(learning_rate = search.best_params_[\"learning_rate\"],\n                           n_estimators  = search.best_params_[\"n_estimators\"],\n                           max_depth     = search.best_params_[\"max_depth\"],\n                           eval_metric='rmsle')\nregressor.fit(X_train9, y_train9)\n\npredictions = regressor.predict(X_test9)\n\nprint('mae is ' , mean_absolute_error(y_test9, predictions))\n\nprint( 'std ground truth is ', y_test9.std())","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:29:09.594582Z","iopub.execute_input":"2022-11-09T12:29:09.595107Z","iopub.status.idle":"2022-11-09T12:37:06.801138Z","shell.execute_reply.started":"2022-11-09T12:29:09.595060Z","shell.execute_reply":"2022-11-09T12:37:06.799794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The model predicts the ordinal feature **position_at_finish** for each horse. For each race, we can rank these predictions to get a model predicted final ranking. For example, if for a race we have the following predictions: {'a': 2, 'b': 2.2, 'c':3.9, 'd': 6.7} where a, b, c and d are program numbers and 2, 2.2, 3.9 and 6.7 are predicted positions at finish, we can rework these numbers to the following ones:  'a': 1, 'b': 2, 'c':3, 'd': 4}","metadata":{}},{"cell_type":"code","source":"preds = pd.Series(predictions, index = y_test9.index, name='prediction')\npositions_preds = pd.merge(y_test9,preds,left_index=True, right_index=True)\npositions_preds.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:06.803285Z","iopub.execute_input":"2022-11-09T12:37:06.804194Z","iopub.status.idle":"2022-11-09T12:37:06.819816Z","shell.execute_reply.started":"2022-11-09T12:37:06.804144Z","shell.execute_reply":"2022-11-09T12:37:06.818085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test9.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:06.821776Z","iopub.execute_input":"2022-11-09T12:37:06.822215Z","iopub.status.idle":"2022-11-09T12:37:06.831965Z","shell.execute_reply.started":"2022-11-09T12:37:06.822180Z","shell.execute_reply":"2022-11-09T12:37:06.830669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test9 = pd.merge(test_df, positions_preds, left_index=True,right_index=True)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:06.833807Z","iopub.execute_input":"2022-11-09T12:37:06.834151Z","iopub.status.idle":"2022-11-09T12:37:06.850018Z","shell.execute_reply.started":"2022-11-09T12:37:06.834120Z","shell.execute_reply":"2022-11-09T12:37:06.848360Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test9.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:06.851889Z","iopub.execute_input":"2022-11-09T12:37:06.852251Z","iopub.status.idle":"2022-11-09T12:37:06.897231Z","shell.execute_reply.started":"2022-11-09T12:37:06.852219Z","shell.execute_reply":"2022-11-09T12:37:06.896046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This code transforms the predicted position into a prediction for each race, as described above, not for each horse. ","metadata":{}},{"cell_type":"code","source":"df_list = []\nfor element, df in test9.groupby(['track_id','race_date','race_number']):\n    df = df.sort_values(by = 'win_prob', ascending=False).reset_index()\n    df['pred_prob'] = df.index + 1\n    df = df.sort_values(by = 'prediction').reset_index()\n    df['pred_mod'] = df.index + 1\n    df_list.append(df)\ndf_all = pd.concat(df_list)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:06.899840Z","iopub.execute_input":"2022-11-09T12:37:06.900346Z","iopub.status.idle":"2022-11-09T12:37:07.817570Z","shell.execute_reply.started":"2022-11-09T12:37:06.900298Z","shell.execute_reply":"2022-11-09T12:37:07.816287Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_prob = accuracy_score(df_all['position_at_finish_x'], df_all['pred_prob'])\nprint(f'The accuracy for predictions derived from the  \\\nprobability of winning alone as calculated from the pre-race odds is {acc_prob: 0.3f}'  )","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:07.819576Z","iopub.execute_input":"2022-11-09T12:37:07.820009Z","iopub.status.idle":"2022-11-09T12:37:07.828035Z","shell.execute_reply.started":"2022-11-09T12:37:07.819973Z","shell.execute_reply":"2022-11-09T12:37:07.826562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_xgb = accuracy_score(df_all['position_at_finish_x'], df_all['pred_mod'])\nprint(f'The accuracy for predictions derived from the  \\\nXGBoost model using only three speed variation features BEFORE the finish, up to the end of the turn, is {acc_xgb: 0.3f}')","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:07.830084Z","iopub.execute_input":"2022-11-09T12:37:07.830687Z","iopub.status.idle":"2022-11-09T12:37:07.844886Z","shell.execute_reply.started":"2022-11-09T12:37:07.830627Z","shell.execute_reply":"2022-11-09T12:37:07.843871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"increase =accuracy_score(df_all['position_at_finish_x'], df_all['pred_mod'])/accuracy_score(df_all['position_at_finish_x'], df_all['pred_prob'])-1\nprint(f'The percentage increase in accuracy of the simple, XGB boost model with features BEFORE the end, up to the turn,  vs odds is {increase: .1%}.' )","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:07.846662Z","iopub.execute_input":"2022-11-09T12:37:07.847407Z","iopub.status.idle":"2022-11-09T12:37:07.859050Z","shell.execute_reply.started":"2022-11-09T12:37:07.847289Z","shell.execute_reply":"2022-11-09T12:37:07.857554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see from the results above, using an XGB classifier on a restricted number of features, only those that describe the race up to the end of the turn, there is an increase in accuracy with respect to the prediction based on the odds alone. ","metadata":{}},{"cell_type":"markdown","source":"### The predictions based on all the features we defined have an accuracy around 84% ","metadata":{}},{"cell_type":"code","source":"X_train9 = train_df[['weight_carried', 'position_at_finish',  'win_prob',\n       'delta_d',  'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts', 'delta_v10', 'delta_v21', 'delta_v32',\n       'delta_d_race', 'draft_v3']]\nX_train9.index.name = 'index'\n\n\ny_train9 = train_df['position_at_finish']\ny_train9.index.name = 'index'\n\ny_test9 = test_df[['position_at_finish']]\ny_test9.index.name ='index'\n\nX_test9 = test_df[['weight_carried', 'position_at_finish',  'win_prob',\n       'delta_d',  'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts', 'delta_v10', 'delta_v21', 'delta_v32',\n       'delta_d_race', 'draft_v3']]\nX_test9.index.name = 'index'\n\n\n\nregressor=xgb.XGBRegressor(eval_metric='rmsle')\nparam_grid = {\"max_depth\":    [4, 5, 6],\n              \"n_estimators\": [500, 600, 700],\n              \"learning_rate\": [0.01, 0.015, 0.02]}\nsearch = GridSearchCV(regressor, param_grid, cv=5).fit(X_train5, y_train5)\n\nprint(\"The best hyperparameters are \",search.best_params_)\nregressor=xgb.XGBRegressor(learning_rate = search.best_params_[\"learning_rate\"],\n                           n_estimators  = search.best_params_[\"n_estimators\"],\n                           max_depth     = search.best_params_[\"max_depth\"],\n                           eval_metric='rmsle')\nregressor.fit(X_train9, y_train9)\n\npredictions = regressor.predict(X_test9)\n\nprint('mae is ' , mean_absolute_error(y_test9, predictions))\n\nprint( 'std ground truth is ', y_test9.std())","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:37:07.860652Z","iopub.execute_input":"2022-11-09T12:37:07.861049Z","iopub.status.idle":"2022-11-09T12:45:07.340950Z","shell.execute_reply.started":"2022-11-09T12:37:07.861014Z","shell.execute_reply":"2022-11-09T12:45:07.339701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds = pd.Series(predictions, index = y_test9.index, name='prediction')\npositions_preds = pd.merge(y_test9,preds,left_index=True, right_index=True)\ntest9 = pd.merge(test_df, positions_preds, left_index=True,right_index=True)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:07.342547Z","iopub.execute_input":"2022-11-09T12:45:07.343630Z","iopub.status.idle":"2022-11-09T12:45:07.354749Z","shell.execute_reply.started":"2022-11-09T12:45:07.343561Z","shell.execute_reply":"2022-11-09T12:45:07.353421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We convert again the \"horse predictions\" to \"race predictions\" by ranking each horse in a race by their **predicted position_at_finish** value.","metadata":{}},{"cell_type":"code","source":"df_list = []\nfor element, df in test9.groupby(['track_id','race_date','race_number']):\n    df = df.sort_values(by = 'win_prob', ascending=False).reset_index()\n    df['pred_prob'] = df.index + 1\n    df = df.sort_values(by = 'prediction').reset_index()\n    df['pred_mod'] = df.index + 1\n    df_list.append(df)\ndf_all = pd.concat(df_list)\ndf_all.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:07.356120Z","iopub.execute_input":"2022-11-09T12:45:07.356489Z","iopub.status.idle":"2022-11-09T12:45:07.765173Z","shell.execute_reply.started":"2022-11-09T12:45:07.356457Z","shell.execute_reply":"2022-11-09T12:45:07.763902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_prob = accuracy_score(df_all['position_at_finish_x'], df_all['pred_prob'])\nprint(f'The accuracy for predictions derived from the  \\\nprobability of winning alone as calculated from the pre-race odds is {acc_prob: 0.3f}'  )","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:07.767080Z","iopub.execute_input":"2022-11-09T12:45:07.767435Z","iopub.status.idle":"2022-11-09T12:45:07.775011Z","shell.execute_reply.started":"2022-11-09T12:45:07.767404Z","shell.execute_reply":"2022-11-09T12:45:07.773593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_xgb = accuracy_score(df_all['position_at_finish_x'], df_all['pred_mod'])\nprint(f'The accuracy for predictions derived from the  \\\nXGBoost model using speed variation features is {acc_xgb: 0.3f}')","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:07.777121Z","iopub.execute_input":"2022-11-09T12:45:07.777632Z","iopub.status.idle":"2022-11-09T12:45:07.796140Z","shell.execute_reply.started":"2022-11-09T12:45:07.777561Z","shell.execute_reply":"2022-11-09T12:45:07.794842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"increase =accuracy_score(df_all['position_at_finish_x'], df_all['pred_mod'])/accuracy_score(df_all['position_at_finish_x'], df_all['pred_prob']) -1\nprint(f'The percentage increase in accuracy of XGB vs odds is {increase: .1%}.' )","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:07.797669Z","iopub.execute_input":"2022-11-09T12:45:07.798141Z","iopub.status.idle":"2022-11-09T12:45:07.810003Z","shell.execute_reply.started":"2022-11-09T12:45:07.798103Z","shell.execute_reply":"2022-11-09T12:45:07.808522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"By using all the features that I defined, the predictions are much better than the predictions based on odds alone. This result confirms that the SHAP values that I got are relevant and we can use those SHAP values to determine which features are more important. ","metadata":{}},{"cell_type":"markdown","source":"# Conclusion ","metadata":{}},{"cell_type":"markdown","source":"The results presented here indicate that the velocity profile along the track is more important than the total path distance (**delta_d**) or the drafting contribution (**draft_v3**). Strategic riding, such as to be close to the rails or to minimize draft **is not worth it**. A better approach is to closely follow an optimal velocity profile for the track and the horse and to minimize velocity fluctuations as much as possible. \nTrack quality and type strongly influences the average horse speed and any analysis on speed or position at finish has to be corrected for the track influence. ","metadata":{}},{"cell_type":"markdown","source":"# Future research","metadata":{}},{"cell_type":"markdown","source":"An interesting future research avenue is to look for better features extracted from the speed variation vs distance along the track. We can also aim to experimentally find an \"ideal speed profile\" along each track and then take the deviation from that profile as another extracted feature. \n\nJust as an indication, I plot the cummulative speed profiles for horses placing 1 and 2 and then for horses placing 5 and above on the same venue, for the same distance. One can see even with the \"naked eye\" that the horses placing above 4, 5 and more have a consistently downward speed profile after the first stretch, as opposed to the horses placing 1 and 2, who can get back to some higher speeds down the race. \n\nFrom the plots below we can think of extracting an \"ideal speed profile\" or \"loser speed profile\" and then use different measures such as mean square distance to find out how close the horse's race is to these profiles. As an observation, because the speed slopes between various points along the track are visibly playing a role, our definitions for extracted features are a good first start. ","metadata":{}},{"cell_type":"code","source":"def make_df(df): # xc2 is the stretch length\n    df_list = []\n    listbad = []\n    for element, frame in tqdm(df.groupby(['track_id','race_date', 'race_number', 'program_number', \\\n        'course_type', 'track_condition', 'run_up_distance', 'race_type', 'purse', 'post_time', 'weight_carried', \\\n        'jockey', 'odds', 'position_at_finish', 'distance', 'win_prob', \\\n                                          \n                                          ])):\n        key = element[0].lower() + '_' + frame.course_type.unique()[0]\n        fp = dict_f[key]\n        xc2 = dict_t[key][1]\n        l = dict_t[key][2]*math.pi\n        \n        try:\n            frame, state = add_vs(frame)\n\n            # here i check wether the horse reached up to the finish line\n            end = frame.loc[frame.pp.transform(lambda z: (z >= fp) & (z.shift(1) < fp))]\n            if len(end) > 0:\n                df_race = frame.query('pp < @fp').sort_values(by=['trakus_index']).reset_index()    \n                df_list.append(df_race)\n        except Exception:\n            listbad.append(element)    \n            \n    return pd.concat(df_list)  ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:45:07.811510Z","iopub.execute_input":"2022-11-09T12:45:07.812481Z","iopub.status.idle":"2022-11-09T12:45:07.824963Z","shell.execute_reply.started":"2022-11-09T12:45:07.812431Z","shell.execute_reply":"2022-11-09T12:45:07.823541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's plot only the sprints at Belmont with distance of 6 furlongs","metadata":{}},{"cell_type":"code","source":"bel_sprint_d = sprints_pp.query('distance_id < 710 & course_type == \"D\" & track_id == \"BEL\"')\nbel_600_12 = bel_sprint_d.query('distance_id == 600 & position_at_finish in [1,2] ')\nbel_600_5 = bel_sprint_d.query('distance_id == 600 & position_at_finish > 4 ')\ndf_12 = make_df(bel_600_12)\ndf_5 = make_df(bel_600_5)","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T12:45:07.827754Z","iopub.execute_input":"2022-11-09T12:45:07.829124Z","iopub.status.idle":"2022-11-09T12:45:21.624680Z","shell.execute_reply.started":"2022-11-09T12:45:07.829072Z","shell.execute_reply":"2022-11-09T12:45:21.623599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import mpl_scatter_density # adds projection='scatter_density'\nfrom matplotlib.colors import LinearSegmentedColormap\n\n# \"Viridis-like\" colormap with white background\nwhite_viridis = LinearSegmentedColormap.from_list('white_viridis', [\n    (0, '#ffffff'),\n    (1e-20, '#440053'),\n    (0.2, '#404388'),\n    (0.4, '#2a788e'),\n    (0.6, '#21a784'),\n    (0.8, '#78d151'),\n    (1, '#fde624'),\n], N=256)\n\ndef using_mpl_scatter_density(fig, x, y):\n    ax = fig.add_subplot(1, 1, 1, projection='scatter_density')\n    density = ax.scatter_density(x, y, cmap=white_viridis)\n    fig.colorbar(density, label='Number of points per pixel')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-09T12:45:21.626926Z","iopub.execute_input":"2022-11-09T12:45:21.627311Z","iopub.status.idle":"2022-11-09T12:45:22.142702Z","shell.execute_reply.started":"2022-11-09T12:45:21.627278Z","shell.execute_reply":"2022-11-09T12:45:22.141350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Density plot for speed vs distance along the track for horses ranking 1 and 2. Notice the bounceback around 1200m and the general keeping of the speed .","metadata":{}},{"cell_type":"code","source":"x = df_12['pp'].to_numpy()\ny = df_12['sps'].to_numpy()\nfig = plt.figure()\nusing_mpl_scatter_density(fig, x, y)\nplt.title(\"Density plot for speed vs distance along the track for horses ranking 1 and 2\")\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T13:27:10.088300Z","iopub.execute_input":"2022-11-09T13:27:10.088792Z","iopub.status.idle":"2022-11-09T13:27:10.412131Z","shell.execute_reply.started":"2022-11-09T13:27:10.088757Z","shell.execute_reply":"2022-11-09T13:27:10.410717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Density plot for speed vs distance along the track for horses ranking 5 and above. Notice that there is little bounceback and the speed does not recover after slowing down at the turn. ","metadata":{}},{"cell_type":"code","source":"x = df_5['pp'].to_numpy()\ny = df_5['sps'].to_numpy()\n\nfig = plt.figure()\nusing_mpl_scatter_density(fig, x, y)\nplt.title(\"Density plot for speed vs distance along the track for horses ranking 5 and above\")\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-09T13:26:40.513797Z","iopub.execute_input":"2022-11-09T13:26:40.514226Z","iopub.status.idle":"2022-11-09T13:26:40.879625Z","shell.execute_reply.started":"2022-11-09T13:26:40.514193Z","shell.execute_reply":"2022-11-09T13:26:40.878030Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Annex I\nThis is just a plot that compares the previously used \"simple draft\" (Spence et al. 2012 ) to the more accurate draft used here, \"draft_v3\"","metadata":{}},{"cell_type":"code","source":"final_corr.plot(x='draft_simple', y='draft_v3', kind = 'scatter')","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:22.907788Z","iopub.execute_input":"2022-11-09T12:45:22.908129Z","iopub.status.idle":"2022-11-09T12:45:23.200188Z","shell.execute_reply.started":"2022-11-09T12:45:22.908099Z","shell.execute_reply":"2022-11-09T12:45:23.198701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Annex II\nHere we make a list of speeds with decreasing standard deviations, all of them summing up to the same value. \nThe last list of speeds is at constant speed, a list of constants. \nThis code snippet is just to show that for a power loss dependence of the form $P \\propto v^{3}$ the best racing strategy is \nconstant speed.\n\n","metadata":{}},{"cell_type":"code","source":"\nindices = np.linspace(1,15,15)\nv3s = []\nfor i in indices:\n    speeds  = (2/i)*np.random.randn(1,200)+20\n    speeds = speeds[0]\n    speeds[199] = speeds[199] - (speeds.sum() - 4000)\n    sv3 = speeds**3\n    v3s.append(sv3.sum())\n\nlast_speed = np.ones(200)*20\nsv3 = last_speed**3\nindices = np.append(indices,16)\nv3s.append(sv3.sum())\nplt.plot(indices, np.asarray(v3s))\nplt.xlabel(\"1/stdev the last point has stdev of 0, constant values\")\nplt.ylabel(\"v^3 summed\")\nplt.title(\"smallest value for the sum of the v^3 function is when the v is constant\")\nplt.show()\n\n\n","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:23.201743Z","iopub.execute_input":"2022-11-09T12:45:23.202150Z","iopub.status.idle":"2022-11-09T12:45:23.392035Z","shell.execute_reply.started":"2022-11-09T12:45:23.202115Z","shell.execute_reply":"2022-11-09T12:45:23.390598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Annex III\nHere I look at one race and see wether the drafting position as I defined it corresponds to the definition. For each horse I plotted both its position relative to the center of mass of all horses (to be able to zoom in) and the position at 5 meters towards the direction of movement. If another horse is within a distance less than these 5 meters and at an angle within plus minus 10 degrees with respect to the direction of movement, then the horse that is \"shadowed\" is in a drafing position. The animation I made is mainly to convince the reader that the function is correct. We can guess the direction of movement from the place we are in the race: left to right in the beginning, right to left at the end and down to up in the middle. \n\nBasically each horse has two points, recognizable by the fact that they have the same shape. Roughly, if there is another point in between those two points, the horse represented by these two points is in drafting position. \n","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"Here I just choose a race","metadata":{}},{"cell_type":"code","source":"counter = 0\nfor element, df in tqdm(sprints_pp.groupby(['track_id','race_date', 'race_number'])):\n        df_race = draft_race(df)\n        counter = counter + 1\n        if counter == 15:\n            break","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:23.393782Z","iopub.execute_input":"2022-11-09T12:45:23.394287Z","iopub.status.idle":"2022-11-09T12:45:42.738934Z","shell.execute_reply.started":"2022-11-09T12:45:23.394249Z","shell.execute_reply":"2022-11-09T12:45:42.737404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" The animation method is taken from this notebook:\n https://www.kaggle.com/code/lichtlab/race-detailed-observation-with-animation","metadata":{}},{"cell_type":"code","source":"color_bool = {True: 'red', False: 'blue'}\nsymbol_bool = {True: 'square', False: 'circle'}","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:42.740951Z","iopub.execute_input":"2022-11-09T12:45:42.741874Z","iopub.status.idle":"2022-11-09T12:45:42.748528Z","shell.execute_reply.started":"2022-11-09T12:45:42.741828Z","shell.execute_reply":"2022-11-09T12:45:42.746976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_av = df_race.groupby('trakus_index').mean()\ndf_av.reset_index(inplace=True)\ndf_av=df_av.rename(columns={'xr':'xrcm','yr':'yrcm'})\ndf_av=df_av[['trakus_index','xrcm','yrcm']]\ndf_haha = pd.merge(df_race,df_av,on='trakus_index')\ndf_haha['x'] = df_haha['xr'] - df_haha['xrcm']\ndf_haha['y'] = df_haha['yr'] - df_haha['yrcm']\ndf_list = []\nfor element, df in df_haha.groupby('program_number'):\n    df['x'] = df['x'] + 5*df['vxs']/(np.sqrt(df['vxs']**2+df['vys']**2))\n    df['y'] = df['y'] + 5*df['vys']/(np.sqrt(df['vxs']**2+df['vys']**2))\n    df_list.append(df)\ndf_plus = pd.concat(df_list)\ndf_all = pd.concat([df_haha,df_plus])","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:42.750383Z","iopub.execute_input":"2022-11-09T12:45:42.750856Z","iopub.status.idle":"2022-11-09T12:45:42.926301Z","shell.execute_reply.started":"2022-11-09T12:45:42.750805Z","shell.execute_reply":"2022-11-09T12:45:42.925009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here I plotted the pairs of points for each horse and whether it is in a drafting position (red) or not (blue). The assignment to a drafting position or not is in line with the definition and what we expect. ","metadata":{}},{"cell_type":"code","source":"px.scatter(df_all, x='x', y='y', animation_frame=\"trakus_index\", \n           color='draft',\n           color_discrete_map=color_bool,\n           symbol='program_number',\n           symbol_map=symbol_bool,\n           range_x = [-10, 10],\n           range_y = [-10, 10],\n           title='Last sprint / .25sec sampling', width=800, height=800)","metadata":{"execution":{"iopub.status.busy":"2022-11-09T12:45:42.928235Z","iopub.execute_input":"2022-11-09T12:45:42.928725Z","iopub.status.idle":"2022-11-09T12:45:51.952654Z","shell.execute_reply.started":"2022-11-09T12:45:42.928684Z","shell.execute_reply":"2022-11-09T12:45:51.951461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Bibliography\n\n**Spence et al. 2012**\nSpeed, pacing strategy and aerodynamic drafting in Thoroughbred horse racing\nAndrew J Spence 1, Andrew S Thurman, Michael J Maher, Alan M Wilson\n*Biol Lett. 2012 Aug 23;8(4):678-81*\n\n**Mercier and Aftalion 2020**\nOptimal speed in Thoroughbred horse racing \nQuentin Mercier and Amandine Aftalion\n*PLoS One. 2020; 15(12): e0235024* \n\n**Davies,  Spence and  Wilson 2019**\nExternal mechanical work in the galloping racehorse\nZ. T. Self Davies, A. J. Spence and A. M. Wilson\n*Biol Lett. 2019 Feb 28;15(2):20180709*\n\n**Pfau ,  Witte and M Wilson 2006**\nCentre of mass movement and mechanical energy fluctuation during gallop locomotion in the Thoroughbred racehorse\nThilo Pfau , Thomas H Witte, Alan M Wilson\n*J Exp Biol. 2006 Oct;209(Pt 19):3742-57*\n\n**Dimensions list**\nUS track dimensions list\nhttps://www.drf.com/sites/default/files/US%20Track%20Dimensions%20List%20-%20May%202015.pdf\n\n\n\n\n","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}