{"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'H:\\Documents\\Kaggle\\big-data-derby-2022'):\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","execution":{"iopub.status.busy":"2023-06-07T11:59:33.388243Z","iopub.execute_input":"2023-06-07T11:59:33.38871Z","iopub.status.idle":"2023-06-07T11:59:33.398462Z","shell.execute_reply.started":"2023-06-07T11:59:33.388673Z","shell.execute_reply":"2023-06-07T11:59:33.396806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Summary\nI have analized 8477 individual horse trajectories that took place in 1134 sprint races (race distance less than 7 furlongs). Contrary to previously published data, our analysis indicates that neither the draft (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 instanteneous speed during the race, indicating that the best strategy for a racing horse is to achieve an optimal speed curve along the track [REF]. 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. ","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":"\n\n### 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\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. \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. 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_{kinetic}}{\\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. $\\newline$\nThe aerodynamic drag also follows a $v^{3}$ dependence (Andrew J. Spence et al., Biology Letters, 2012, 8, 678-681): $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 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}$\n\nNow that we have this proportionality relationthip we can make the following observation: to minimize the total 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. \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, 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}$, therefor the total 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. ","metadata":{}},{"cell_type":"markdown","source":"### Percentage effects of aerodynamic draft and path efficiency are small.\n\nExisting literature puts the contribution of aerodynamic drag to the total energy expenditure during the race at 17%. (DAvies, Spence, Bilogy Letters). If the aerodynamic drafting is reduces by a factor of say, 20%, due to \"strategic riding\", the energy saved would be around 3.4%. For the estimation of the aerodynamic 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 for the drag is an upper limit and the actual value is lower. Experiments in wind tunnels are needed to measure exactly the aerodynamic 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 bing 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 drafting 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 maximizing aerodynamic drag are worth it. ","metadata":{}},{"cell_type":"markdown","source":"Let's 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":"We begin by transforming the latitude and longitude coordinates into $x$ and $y$ coordinates. The $x$ and $y$ coordinates will then be rotated and and translated so that the chute is at negative $x$ values, the chute and the first stretch are paralel to the $x$ axis and the center of the leftmost turn (the one above the chute) is at position $ x= 0$ and $y=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. The actual projected position is the position that counts in terms of advancing along the track during the race. \n\nUsing google maps I have manually got the latitude and longitude points of the two lines (chute + first stretch & final stretch) and the two turns. These will be transformed in x and y variables and then fit to a line or circle respectively. Once we get the lines and circles we can rotate and translate the coordinates axes to get the optimal axes for determining the projected position (pp) for each point in the track. ","metadata":{}},{"cell_type":"code","source":"!pip install pymap3d\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\nfrom tqdm import tqdm\nimport plotly.express as px\n#plt.style.use('ggplot')\n#import mpl_scatter_density # adds projection='scatter_density'\nfrom matplotlib.colors import LinearSegmentedColormap","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:33.400972Z","iopub.execute_input":"2023-06-07T11:59:33.401753Z","iopub.status.idle":"2023-06-07T11:59:46.657476Z","shell.execute_reply.started":"2023-06-07T11:59:33.401715Z","shell.execute_reply":"2023-06-07T11:59:46.656143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os.path","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.659352Z","iopub.execute_input":"2023-06-07T11:59:46.659698Z","iopub.status.idle":"2023-06-07T11:59:46.665618Z","shell.execute_reply.started":"2023-06-07T11:59:46.659663Z","shell.execute_reply":"2023-06-07T11:59:46.664447Z"},"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","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.669458Z","iopub.execute_input":"2023-06-07T11:59:46.669854Z","iopub.status.idle":"2023-06-07T11:59:46.679752Z","shell.execute_reply.started":"2023-06-07T11:59:46.66982Z","shell.execute_reply":"2023-06-07T11:59:46.678803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sklearn","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.681303Z","iopub.execute_input":"2023-06-07T11:59:46.681626Z","iopub.status.idle":"2023-06-07T11:59:46.693403Z","shell.execute_reply.started":"2023-06-07T11:59:46.681595Z","shell.execute_reply":"2023-06-07T11:59:46.692479Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings # Supress warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.695084Z","iopub.execute_input":"2023-06-07T11:59:46.695773Z","iopub.status.idle":"2023-06-07T11:59:46.706279Z","shell.execute_reply.started":"2023-06-07T11:59:46.695734Z","shell.execute_reply":"2023-06-07T11:59:46.70527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BASE_DIR=r'/kaggle/input/big-data-derby-2022'","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.707965Z","iopub.execute_input":"2023-06-07T11:59:46.708407Z","iopub.status.idle":"2023-06-07T11:59:46.718928Z","shell.execute_reply.started":"2023-06-07T11:59:46.708374Z","shell.execute_reply":"2023-06-07T11:59:46.717766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xp1 = os.path.join(\"data\",\"cities\",\"geo.txt\")","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.720655Z","iopub.execute_input":"2023-06-07T11:59:46.721403Z","iopub.status.idle":"2023-06-07T11:59:46.733356Z","shell.execute_reply.started":"2023-06-07T11:59:46.72137Z","shell.execute_reply":"2023-06-07T11:59:46.732305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### How to get the exact position of the racetrack\nEach track on the racetrack is composed of two line segments and two semicircles. We can fit them to points collected by hand. I have collected by hand using google maps a collection of points for each segment and each semicirfle 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\". ","metadata":{}},{"cell_type":"markdown","source":"Each racetrack has three tracks, the D, the turf T and the inner turf I. The names for the Aqueduct are different and we will take that aspect into consideration. The satellite image for the Aqueduct inner turf are not good in the sense that I cannot see the rails exactly so I left out the Aqueduct inner turf races. 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":"code","source":"# ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.735167Z","iopub.execute_input":"2023-06-07T11:59:46.735588Z","iopub.status.idle":"2023-06-07T11:59:46.745949Z","shell.execute_reply.started":"2023-06-07T11:59:46.735556Z","shell.execute_reply":"2023-06-07T11:59:46.744917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# points to fit the venues\nall_points = pd.read_csv(r'/kaggle/input/points-for-fitting-tracks/all_points.csv')\nall_points.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.750775Z","iopub.execute_input":"2023-06-07T11:59:46.751534Z","iopub.status.idle":"2023-06-07T11:59:46.773625Z","shell.execute_reply.started":"2023-06-07T11:59:46.751499Z","shell.execute_reply":"2023-06-07T11:59:46.77229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_points.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.775472Z","iopub.execute_input":"2023-06-07T11:59:46.77642Z","iopub.status.idle":"2023-06-07T11:59:46.792716Z","shell.execute_reply.started":"2023-06-07T11:59:46.776383Z","shell.execute_reply":"2023-06-07T11:59:46.791642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The finish line as identified via google maps\naqu_f = {'latitude' : [40.671707525293485], \"longitude\": [-73.83180709491133], 'course_type' : ['finish'], 'track_id' : [\"AQU\"], 'distance_id' : [600]}\ndf_l = pd.DataFrame(aqu_f)\ndf_l","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.796505Z","iopub.execute_input":"2023-06-07T11:59:46.797307Z","iopub.status.idle":"2023-06-07T11:59:46.814736Z","shell.execute_reply.started":"2023-06-07T11:59:46.797268Z","shell.execute_reply":"2023-06-07T11:59:46.813467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"files_dict = {}\nfor dirname, _, filenames in os.walk(BASE_DIR):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n        files_dict[filename] = pd.read_csv(os.path.join(dirname, filename))","metadata":{"execution":{"iopub.status.busy":"2023-06-07T11:59:46.816809Z","iopub.execute_input":"2023-06-07T11:59:46.817922Z","iopub.status.idle":"2023-06-07T12:00:09.854826Z","shell.execute_reply.started":"2023-06-07T11:59:46.817878Z","shell.execute_reply":"2023-06-07T12:00:09.853638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"files_list=[]\nfor key in files_dict.keys():\n    print(key)\n    files_list.append(key)\nprint(files_list)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.856481Z","iopub.execute_input":"2023-06-07T12:00:09.856938Z","iopub.status.idle":"2023-06-07T12:00:09.865274Z","shell.execute_reply.started":"2023-06-07T12:00:09.856894Z","shell.execute_reply":"2023-06-07T12:00:09.864055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race = files_dict['nyra_race_table.csv']","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.867518Z","iopub.execute_input":"2023-06-07T12:00:09.867958Z","iopub.status.idle":"2023-06-07T12:00:09.884077Z","shell.execute_reply.started":"2023-06-07T12:00:09.867918Z","shell.execute_reply":"2023-06-07T12:00:09.882318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.885867Z","iopub.execute_input":"2023-06-07T12:00:09.886297Z","iopub.status.idle":"2023-06-07T12:00:09.905006Z","shell.execute_reply.started":"2023-06-07T12:00:09.886263Z","shell.execute_reply":"2023-06-07T12:00:09.903689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.906832Z","iopub.execute_input":"2023-06-07T12:00:09.907254Z","iopub.status.idle":"2023-06-07T12:00:09.931916Z","shell.execute_reply.started":"2023-06-07T12:00:09.907212Z","shell.execute_reply":"2023-06-07T12:00:09.930371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race[['year', 'month', 'day']] = df_race['race_date'].str.split('-', expand=True)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.933746Z","iopub.execute_input":"2023-06-07T12:00:09.934131Z","iopub.status.idle":"2023-06-07T12:00:09.956674Z","shell.execute_reply.started":"2023-06-07T12:00:09.934073Z","shell.execute_reply":"2023-06-07T12:00:09.95513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race['year'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.95884Z","iopub.execute_input":"2023-06-07T12:00:09.959267Z","iopub.status.idle":"2023-06-07T12:00:09.972345Z","shell.execute_reply.started":"2023-06-07T12:00:09.959222Z","shell.execute_reply":"2023-06-07T12:00:09.971146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race['furlongs'] = df_race['distance_id']/100","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.973756Z","iopub.execute_input":"2023-06-07T12:00:09.974106Z","iopub.status.idle":"2023-06-07T12:00:09.98615Z","shell.execute_reply.started":"2023-06-07T12:00:09.974077Z","shell.execute_reply":"2023-06-07T12:00:09.984905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"long_races = df_race.query('furlongs > 7.05')\nlong_races.shape\n                           ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:09.987914Z","iopub.execute_input":"2023-06-07T12:00:09.988546Z","iopub.status.idle":"2023-06-07T12:00:10.010394Z","shell.execute_reply.started":"2023-06-07T12:00:09.988498Z","shell.execute_reply":"2023-06-07T12:00:10.008886Z"},"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. In our dataset the short races are the most numerous more than half the total number of races and I will focus on them for now. ","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":"2023-06-07T12:00:10.012432Z","iopub.execute_input":"2023-06-07T12:00:10.012789Z","iopub.status.idle":"2023-06-07T12:00:10.468479Z","shell.execute_reply.started":"2023-06-07T12:00:10.01276Z","shell.execute_reply":"2023-06-07T12:00:10.46727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Following [REF] we are going to analyse sprints in our work, that is races 7 furlongs and less. ","metadata":{}},{"cell_type":"code","source":"df_race['furlongs'].min()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:10.470456Z","iopub.execute_input":"2023-06-07T12:00:10.470958Z","iopub.status.idle":"2023-06-07T12:00:10.480812Z","shell.execute_reply.started":"2023-06-07T12:00:10.470913Z","shell.execute_reply":"2023-06-07T12:00:10.479317Z"},"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)\nprint(f'The percentage of short races is {percentage: .1%}.' )\n                       ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:10.48264Z","iopub.execute_input":"2023-06-07T12:00:10.482993Z","iopub.status.idle":"2023-06-07T12:00:10.502528Z","shell.execute_reply.started":"2023-06-07T12:00:10.482963Z","shell.execute_reply":"2023-06-07T12:00:10.501117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"short_races.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:10.506914Z","iopub.execute_input":"2023-06-07T12:00:10.507322Z","iopub.status.idle":"2023-06-07T12:00:10.520696Z","shell.execute_reply.started":"2023-06-07T12:00:10.507289Z","shell.execute_reply":"2023-06-07T12:00:10.519394Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Here I load the complete dataset","metadata":{}},{"cell_type":"code","source":"df_complete = files_dict['nyra_2019_complete.csv'].copy()\ndf_complete.columns","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:10.522305Z","iopub.execute_input":"2023-06-07T12:00:10.522767Z","iopub.status.idle":"2023-06-07T12:00:11.44076Z","shell.execute_reply.started":"2023-06-07T12:00:10.522732Z","shell.execute_reply":"2023-06-07T12:00:11.439467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_complete.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:11.442755Z","iopub.execute_input":"2023-06-07T12:00:11.443239Z","iopub.status.idle":"2023-06-07T12:00:11.450282Z","shell.execute_reply.started":"2023-06-07T12:00:11.443173Z","shell.execute_reply":"2023-06-07T12:00:11.449376Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:11.461637Z","iopub.execute_input":"2023-06-07T12:00:11.462111Z","iopub.status.idle":"2023-06-07T12:00:12.309796Z","shell.execute_reply.started":"2023-06-07T12:00:11.46207Z","shell.execute_reply":"2023-06-07T12:00:12.308324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_complete.course_type.unique()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:12.312139Z","iopub.execute_input":"2023-06-07T12:00:12.3127Z","iopub.status.idle":"2023-06-07T12:00:12.740057Z","shell.execute_reply.started":"2023-06-07T12:00:12.312639Z","shell.execute_reply":"2023-06-07T12:00:12.738834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# here I add the points that I will use to fit the racetracks \ndf_complete = pd.concat([df_complete, all_points])","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:12.741784Z","iopub.execute_input":"2023-06-07T12:00:12.74216Z","iopub.status.idle":"2023-06-07T12:00:13.43146Z","shell.execute_reply.started":"2023-06-07T12:00:12.742126Z","shell.execute_reply":"2023-06-07T12:00:13.430276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_complete.columns","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:13.433154Z","iopub.execute_input":"2023-06-07T12:00:13.433651Z","iopub.status.idle":"2023-06-07T12:00:13.442052Z","shell.execute_reply.started":"2023-06-07T12:00:13.433599Z","shell.execute_reply":"2023-06-07T12:00:13.440801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_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":"2023-06-07T12:00:13.444041Z","iopub.execute_input":"2023-06-07T12:00:13.444517Z","iopub.status.idle":"2023-06-07T12:00:17.221216Z","shell.execute_reply.started":"2023-06-07T12:00:13.444475Z","shell.execute_reply":"2023-06-07T12:00:17.219694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_complete['distance'] = (df_complete['distance_id']/100)*201.168 + df_complete['run_up_distance']*0.3048","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:17.223022Z","iopub.execute_input":"2023-06-07T12:00:17.223725Z","iopub.status.idle":"2023-06-07T12:00:17.288001Z","shell.execute_reply.started":"2023-06-07T12:00:17.223685Z","shell.execute_reply":"2023-06-07T12:00:17.28673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_complete['odds'] = df_complete['odds']/100","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:17.290809Z","iopub.execute_input":"2023-06-07T12:00:17.292156Z","iopub.status.idle":"2023-06-07T12:00:17.358533Z","shell.execute_reply.started":"2023-06-07T12:00:17.2921Z","shell.execute_reply":"2023-06-07T12:00:17.357412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_complete['win_prob'] =  1/(df_complete['odds'] + 1)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:17.361073Z","iopub.execute_input":"2023-06-07T12:00:17.36189Z","iopub.status.idle":"2023-06-07T12:00:17.428709Z","shell.execute_reply.started":"2023-06-07T12:00:17.36183Z","shell.execute_reply":"2023-06-07T12:00:17.427671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Following Joe Appelbaum's advice in teh discussions I will look only at course_type D, O, I, T and leave M aside. ","metadata":{}},{"cell_type":"code","source":"aqu = df_complete.query('track_id == \"AQU\"')","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:17.430005Z","iopub.execute_input":"2023-06-07T12:00:17.431029Z","iopub.status.idle":"2023-06-07T12:00:19.889499Z","shell.execute_reply.started":"2023-06-07T12:00:17.430989Z","shell.execute_reply":"2023-06-07T12:00:19.888261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.title('Course type pie')\nplt.pie(df_complete['course_type'].value_counts(), labels=df_complete['course_type'].unique())\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:19.891098Z","iopub.execute_input":"2023-06-07T12:00:19.891546Z","iopub.status.idle":"2023-06-07T12:00:21.260494Z","shell.execute_reply.started":"2023-06-07T12:00:19.891506Z","shell.execute_reply":"2023-06-07T12:00:21.259222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bel = df_complete.query('track_id == \"BEL\"')","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:21.262125Z","iopub.execute_input":"2023-06-07T12:00:21.263472Z","iopub.status.idle":"2023-06-07T12:00:21.875286Z","shell.execute_reply.started":"2023-06-07T12:00:21.263411Z","shell.execute_reply":"2023-06-07T12:00:21.874032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.title('Course type pie BEL')\nplt.pie(bel['course_type'].value_counts(), labels=bel['course_type'].unique())\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:21.877095Z","iopub.execute_input":"2023-06-07T12:00:21.877485Z","iopub.status.idle":"2023-06-07T12:00:22.458817Z","shell.execute_reply.started":"2023-06-07T12:00:21.877453Z","shell.execute_reply":"2023-06-07T12:00:22.457625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sar = df_complete.query('track_id == \"SAR\"')\nplt.pie(sar['course_type'].value_counts(), labels=sar['course_type'].unique())\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:22.460757Z","iopub.execute_input":"2023-06-07T12:00:22.461546Z","iopub.status.idle":"2023-06-07T12:00:23.342087Z","shell.execute_reply.started":"2023-06-07T12:00:22.4615Z","shell.execute_reply":"2023-06-07T12:00:23.340966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### We can get 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> -   transform the (latitude,longitude) coordinates into (x,y,z) coordinates using the pymap3d python module\n> -  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 and getting the transformation matrix from the eigenvectors of the matrix containing the moments of inertia. \n> -  fit the lines corresponding to the two stretches \n> -   rotate the axes of coordinates so that the two stretches are horizontal\n> -   translate the axes of coordinates so that the center of the leftmost semicircle is at (0,0). \n>  These coordinates are designates with xr, yr. \n","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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:23.344409Z","iopub.execute_input":"2023-06-07T12:00:23.345206Z","iopub.status.idle":"2023-06-07T12:00:23.362666Z","shell.execute_reply.started":"2023-06-07T12:00:23.345144Z","shell.execute_reply":"2023-06-07T12:00:23.361079Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:23.364831Z","iopub.execute_input":"2023-06-07T12:00:23.365359Z","iopub.status.idle":"2023-06-07T12:00:23.38953Z","shell.execute_reply.started":"2023-06-07T12:00:23.365312Z","shell.execute_reply":"2023-06-07T12:00:23.387965Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:23.391271Z","iopub.execute_input":"2023-06-07T12:00:23.39226Z","iopub.status.idle":"2023-06-07T12:00:23.407737Z","shell.execute_reply.started":"2023-06-07T12:00:23.392209Z","shell.execute_reply":"2023-06-07T12:00:23.406762Z"},"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('x (meters)')\n    plt.ylabel('y (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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:23.408926Z","iopub.execute_input":"2023-06-07T12:00:23.409819Z","iopub.status.idle":"2023-06-07T12:00:23.436218Z","shell.execute_reply.started":"2023-06-07T12:00:23.409785Z","shell.execute_reply":"2023-06-07T12:00:23.434585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#The adding of xr, yr coordinates, the good rotated coordinates have to be done for each venue\n# at a time\naqu = df_complete.query('track_id == \"AQU\" ')\naqu = add_xy(aqu)\nprint('shape aqu sprint ', aqu.shape)\nbel = df_complete.query('track_id == \"BEL\" ')\nbel = add_xy(bel)\nprint('shape bel sprint ', bel.shape)\nsar = df_complete.query('track_id == \"SAR\" ')\nsar = add_xy_sar(sar)\nprint('shape sar sprint ', sar.shape)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:23.437533Z","iopub.execute_input":"2023-06-07T12:00:23.437857Z","iopub.status.idle":"2023-06-07T12:00:29.430564Z","shell.execute_reply.started":"2023-06-07T12:00:23.43783Z","shell.execute_reply":"2023-06-07T12:00:29.429109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# in this dictionary I put all the dataframes of the rotated and translated coordinates\n# together with info on fits:  radius of the turn, length of the stretch and translation distance\ndict_t={}","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:29.431842Z","iopub.execute_input":"2023-06-07T12:00:29.43221Z","iopub.status.idle":"2023-06-07T12:00:29.438751Z","shell.execute_reply.started":"2023-06-07T12:00:29.432157Z","shell.execute_reply":"2023-06-07T12:00:29.437119Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:29.440012Z","iopub.execute_input":"2023-06-07T12:00:29.441324Z","iopub.status.idle":"2023-06-07T12:00:36.847883Z","shell.execute_reply.started":"2023-06-07T12:00:29.441288Z","shell.execute_reply":"2023-06-07T12:00:36.846551Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:36.849641Z","iopub.execute_input":"2023-06-07T12:00:36.850065Z","iopub.status.idle":"2023-06-07T12:00:41.427605Z","shell.execute_reply.started":"2023-06-07T12:00:36.850027Z","shell.execute_reply":"2023-06-07T12:00:41.426445Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:41.429251Z","iopub.execute_input":"2023-06-07T12:00:41.429658Z","iopub.status.idle":"2023-06-07T12:00:48.993817Z","shell.execute_reply.started":"2023-06-07T12:00:41.429624Z","shell.execute_reply":"2023-06-07T12:00:48.992848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_t.keys()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:48.995415Z","iopub.execute_input":"2023-06-07T12:00:48.99605Z","iopub.status.idle":"2023-06-07T12:00:49.002768Z","shell.execute_reply.started":"2023-06-07T12:00:48.996016Z","shell.execute_reply":"2023-06-07T12:00:49.001298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# (1 stretch, 2 radius, 3 translation_x)\ndict_t['bel_T'][1:]","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:49.004309Z","iopub.execute_input":"2023-06-07T12:00:49.005067Z","iopub.status.idle":"2023-06-07T12:00:49.02129Z","shell.execute_reply.started":"2023-06-07T12:00:49.005022Z","shell.execute_reply":"2023-06-07T12:00:49.019895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_t['aqu_D'][1:]","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:49.023148Z","iopub.execute_input":"2023-06-07T12:00:49.02397Z","iopub.status.idle":"2023-06-07T12:00:49.037821Z","shell.execute_reply.started":"2023-06-07T12:00:49.023933Z","shell.execute_reply":"2023-06-07T12:00:49.036359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### The track dimensions obtained from the fits are similar to the ones indicated in the venue descriptions","metadata":{}},{"cell_type":"markdown","source":" The 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":{}},{"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. ","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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:49.039809Z","iopub.execute_input":"2023-06-07T12:00:49.040226Z","iopub.status.idle":"2023-06-07T12:00:49.055524Z","shell.execute_reply.started":"2023-06-07T12:00:49.040167Z","shell.execute_reply":"2023-06-07T12:00:49.054302Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:49.05715Z","iopub.execute_input":"2023-06-07T12:00:49.057529Z","iopub.status.idle":"2023-06-07T12:00:56.321643Z","shell.execute_reply.started":"2023-06-07T12:00:49.057498Z","shell.execute_reply":"2023-06-07T12:00:56.320303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sprints_pp = sprints_pp.drop(columns=['latitude','longitude','x','y','z'])","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:56.322872Z","iopub.execute_input":"2023-06-07T12:00:56.323286Z","iopub.status.idle":"2023-06-07T12:00:56.725961Z","shell.execute_reply.started":"2023-06-07T12:00:56.323248Z","shell.execute_reply":"2023-06-07T12:00:56.724842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sprints_pp.columns","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:56.727731Z","iopub.execute_input":"2023-06-07T12:00:56.728071Z","iopub.status.idle":"2023-06-07T12:00:56.735889Z","shell.execute_reply.started":"2023-06-07T12:00:56.728042Z","shell.execute_reply":"2023-06-07T12:00:56.734678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2023-06-07T12:00:56.737503Z","iopub.execute_input":"2023-06-07T12:00:56.737934Z","iopub.status.idle":"2023-06-07T12:00:56.754435Z","shell.execute_reply.started":"2023-06-07T12:00:56.737893Z","shell.execute_reply":"2023-06-07T12:00:56.75296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sprints_pp.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:56.756067Z","iopub.execute_input":"2023-06-07T12:00:56.756565Z","iopub.status.idle":"2023-06-07T12:00:56.779327Z","shell.execute_reply.started":"2023-06-07T12:00:56.756532Z","shell.execute_reply":"2023-06-07T12:00:56.77802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# I just check that all the chutes are having a shape that this model understands. \nsprints_pp.plot(x='xr', y='yr', kind='scatter' )","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:00:56.780967Z","iopub.execute_input":"2023-06-07T12:00:56.781433Z","iopub.status.idle":"2023-06-07T12:01:06.8318Z","shell.execute_reply.started":"2023-06-07T12:00:56.781386Z","shell.execute_reply":"2023-06-07T12:01:06.830378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_t['bel_D'][3]","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.83336Z","iopub.execute_input":"2023-06-07T12:01:06.833725Z","iopub.status.idle":"2023-06-07T12:01:06.843278Z","shell.execute_reply.started":"2023-06-07T12:01:06.833695Z","shell.execute_reply":"2023-06-07T12:01:06.841867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Aqueduct inner track (outer turf)\n1 Mile Around\n1174.9' Head of Stretch to Finish Line\n332.6' Finish Line to First Turn\n350' Turn Radius\n1132.5 +' Length of Turn","metadata":{}},{"cell_type":"markdown","source":"BELMONT\nMain Track: 1 1/2 Miles Around\n1097' Head of Stretch to Finish Line\n843' Finish Line to First Turn\n643' Turn Radius\n2020' Length of Turn","metadata":{}},{"cell_type":"code","source":"# these are the finishing 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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.845143Z","iopub.execute_input":"2023-06-07T12:01:06.845894Z","iopub.status.idle":"2023-06-07T12:01:06.856155Z","shell.execute_reply.started":"2023-06-07T12:01:06.845848Z","shell.execute_reply":"2023-06-07T12:01:06.854762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_f","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.858005Z","iopub.execute_input":"2023-06-07T12:01:06.858481Z","iopub.status.idle":"2023-06-07T12:01:06.870901Z","shell.execute_reply.started":"2023-06-07T12:01:06.858437Z","shell.execute_reply":"2023-06-07T12:01:06.869289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# let's check the distance around at Belmont Turf\n# Widener Turf: 1 5/16 Miles + 27' Around\n# in meters from the published track dimensions\n2112.264 + 8.2296","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.872382Z","iopub.execute_input":"2023-06-07T12:01:06.872734Z","iopub.status.idle":"2023-06-07T12:01:06.88605Z","shell.execute_reply.started":"2023-06-07T12:01:06.872702Z","shell.execute_reply":"2023-06-07T12:01:06.88447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#From our fits the track distance around at Belmont Turf. A 4 metre deviation, 0.1%\n2*dict_t['bel_T'][1]+2*dict_t['bel_T'][2]*math.pi","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.88789Z","iopub.execute_input":"2023-06-07T12:01:06.88834Z","iopub.status.idle":"2023-06-07T12:01:06.903481Z","shell.execute_reply.started":"2023-06-07T12:01:06.888306Z","shell.execute_reply":"2023-06-07T12:01:06.902018Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.905135Z","iopub.execute_input":"2023-06-07T12:01:06.905636Z","iopub.status.idle":"2023-06-07T12:01:06.930154Z","shell.execute_reply.started":"2023-06-07T12:01:06.905591Z","shell.execute_reply":"2023-06-07T12:01:06.928699Z"},"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] ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.931644Z","iopub.execute_input":"2023-06-07T12:01:06.931993Z","iopub.status.idle":"2023-06-07T12:01:06.949003Z","shell.execute_reply.started":"2023-06-07T12:01:06.931956Z","shell.execute_reply":"2023-06-07T12:01:06.94796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_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] ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.950589Z","iopub.execute_input":"2023-06-07T12:01:06.950948Z","iopub.status.idle":"2023-06-07T12:01:06.96532Z","shell.execute_reply.started":"2023-06-07T12:01:06.950918Z","shell.execute_reply":"2023-06-07T12:01:06.963538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_t['bel_D'][3] ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.967011Z","iopub.execute_input":"2023-06-07T12:01:06.96741Z","iopub.status.idle":"2023-06-07T12:01:06.984622Z","shell.execute_reply.started":"2023-06-07T12:01:06.967375Z","shell.execute_reply":"2023-06-07T12:01:06.983629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Saratoga Turf:\n#Mellon Turf: 1 Mile Around\n#~380' Turn Radius","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:06.986033Z","iopub.execute_input":"2023-06-07T12:01:06.986434Z","iopub.status.idle":"2023-06-07T12:01:07.002795Z","shell.execute_reply.started":"2023-06-07T12:01:06.986399Z","shell.execute_reply":"2023-06-07T12:01:07.001541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"SARATOGA\nMain Track: 1 1/8 Miles Around\n1144' Head of Stretch to Finish Line\n341' Finish Line to First Turn\n~470' Turn Radius\n1485' Length of Turn\nMellon Turf: 1 Mile Around\n~380' Turn Radius\nInner Turf / Steeplechase: 7 Furlongs Around","metadata":{}},{"cell_type":"code","source":"dict_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] ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:07.004393Z","iopub.execute_input":"2023-06-07T12:01:07.004834Z","iopub.status.idle":"2023-06-07T12:01:07.019466Z","shell.execute_reply.started":"2023-06-07T12:01:07.004794Z","shell.execute_reply":"2023-06-07T12:01:07.017993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:07.02133Z","iopub.execute_input":"2023-06-07T12:01:07.021795Z","iopub.status.idle":"2023-06-07T12:01:07.036468Z","shell.execute_reply.started":"2023-06-07T12:01:07.021751Z","shell.execute_reply":"2023-06-07T12:01:07.034975Z"},"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":"2023-06-07T12:01:07.038136Z","iopub.execute_input":"2023-06-07T12:01:07.038633Z","iopub.status.idle":"2023-06-07T12:01:07.052845Z","shell.execute_reply.started":"2023-06-07T12:01:07.038586Z","shell.execute_reply":"2023-06-07T12:01:07.051429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This function adds the actual speed, 'sa' and the projected speed 'sp' to the dataframe. It smoothes the speed obtained from the data by applying a savgol filter as in the PLOS ONE paper. 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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:07.054863Z","iopub.execute_input":"2023-06-07T12:01:07.05532Z","iopub.status.idle":"2023-06-07T12:01:07.072135Z","shell.execute_reply.started":"2023-06-07T12:01:07.055277Z","shell.execute_reply":"2023-06-07T12:01:07.070871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The code below shows, for illustration, the actual speed and the projected speed (lower) for a horse chosen at random. ","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","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:07.073904Z","iopub.execute_input":"2023-06-07T12:01:07.075617Z","iopub.status.idle":"2023-06-07T12:01:08.561306Z","shell.execute_reply.started":"2023-06-07T12:01:07.075575Z","shell.execute_reply":"2023-06-07T12:01:08.559883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"datenow = 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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:08.563092Z","iopub.execute_input":"2023-06-07T12:01:08.564986Z","iopub.status.idle":"2023-06-07T12:01:08.901962Z","shell.execute_reply.started":"2023-06-07T12:01:08.564934Z","shell.execute_reply":"2023-06-07T12:01:08.900456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is a plot of the actual speed and the projected speed. Notice how a the turn the actual speed is higher than the \nprojected speed as expected.","metadata":{}},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (12,12)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:08.903812Z","iopub.execute_input":"2023-06-07T12:01:08.904275Z","iopub.status.idle":"2023-06-07T12:01:08.909551Z","shell.execute_reply.started":"2023-06-07T12:01:08.904239Z","shell.execute_reply":"2023-06-07T12:01:08.908307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:08.91115Z","iopub.execute_input":"2023-06-07T12:01:08.911568Z","iopub.status.idle":"2023-06-07T12:01:09.361417Z","shell.execute_reply.started":"2023-06-07T12:01:08.911536Z","shell.execute_reply":"2023-06-07T12:01:09.360176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## We define a few quantities that we will extract 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 \n#### Energy minimization features\n> - **draft percentage** is the percentage of the aerodynamic drag that is minimized during the race. We followed [REF] in defining the draft percentage 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}$. The draft percentage I used takes into account that being shielded at low sppeds 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 = \\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 behind at least another horse at an angle no larger than $10^{o}$. The formula used here and the formula used in [REF] give the same result at 0 and 1. \n> The draft percentage used in [REF] 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## The question is: which type of features is more important in the outcome of an individual race? The outcomes are projected average speed (pa) and positon at finish.","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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:09.363354Z","iopub.execute_input":"2023-06-07T12:01:09.364075Z","iopub.status.idle":"2023-06-07T12:01:09.38726Z","shell.execute_reply.started":"2023-06-07T12:01:09.364038Z","shell.execute_reply":"2023-06-07T12:01:09.385838Z"},"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    ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:09.389435Z","iopub.execute_input":"2023-06-07T12:01:09.390174Z","iopub.status.idle":"2023-06-07T12:01:09.41068Z","shell.execute_reply.started":"2023-06-07T12:01:09.390086Z","shell.execute_reply":"2023-06-07T12:01:09.409681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"limit = math.cos(3.1415*10/180)\nlimit","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:09.41206Z","iopub.execute_input":"2023-06-07T12:01:09.412786Z","iopub.status.idle":"2023-06-07T12:01:09.434468Z","shell.execute_reply.started":"2023-06-07T12:01:09.41275Z","shell.execute_reply":"2023-06-07T12:01:09.433264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bool(2)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:09.436096Z","iopub.execute_input":"2023-06-07T12:01:09.436895Z","iopub.status.idle":"2023-06-07T12:01:09.45083Z","shell.execute_reply.started":"2023-06-07T12:01:09.43685Z","shell.execute_reply":"2023-06-07T12:01:09.449262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df[['race_number','distance_id']].sum(axis=1).std()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:09.478725Z","iopub.execute_input":"2023-06-07T12:01:09.479211Z","iopub.status.idle":"2023-06-07T12:01:09.502901Z","shell.execute_reply.started":"2023-06-07T12:01:09.479162Z","shell.execute_reply":"2023-06-07T12:01:09.501671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race.columns","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:09.507089Z","iopub.execute_input":"2023-06-07T12:01:09.507757Z","iopub.status.idle":"2023-06-07T12:01:09.514927Z","shell.execute_reply.started":"2023-06-07T12:01:09.507718Z","shell.execute_reply":"2023-06-07T12:01:09.513997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":"2023-06-07T12:01:09.516289Z","iopub.execute_input":"2023-06-07T12:01:09.51685Z","iopub.status.idle":"2023-06-07T12:01:29.479572Z","shell.execute_reply.started":"2023-06-07T12:01:09.516817Z","shell.execute_reply":"2023-06-07T12:01:29.478297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/code/lichtlab/race-detailed-observation-with-animation","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:29.481077Z","iopub.execute_input":"2023-06-07T12:01:29.481994Z","iopub.status.idle":"2023-06-07T12:01:29.486959Z","shell.execute_reply.started":"2023-06-07T12:01:29.481958Z","shell.execute_reply":"2023-06-07T12:01:29.485652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"color_discrete_map = {\n    \"6\": \"red\",\n    \"1\": \"green\",\n    \"11\": \"blue\",\n    \"9\": \"goldenrod\",\n    \"4\": \"magenta\",\n    \"7\": \"orange\",\n    \"10\": \"pink\",\n    \"3\": \"gray\",\n    \"2\": \"gray\",\n    \"8\": \"gray\"\n}\n\nsymbol_discrete_map = {\n    \"6\": \"star\",\n    \"1\": \"square\",\n    \"11\": \"x\",\n    \"9\": \"circle\",\n    \"4\": \"cross\",\n    \"7\": \"diamond\",\n    \"10\": \"triangle-up\",\n    \"3\": \"circle-x\",\n    \"2\": \"circle-x\",\n    \"8\": \"circle-x\"    \n}\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:29.488751Z","iopub.execute_input":"2023-06-07T12:01:29.489091Z","iopub.status.idle":"2023-06-07T12:01:29.504348Z","shell.execute_reply.started":"2023-06-07T12:01:29.489062Z","shell.execute_reply":"2023-06-07T12:01:29.503043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"color_bool = {True: 'red', False: 'blue'}\nsymbol_bool = {True: 'square', False: 'circle'}","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:29.506358Z","iopub.execute_input":"2023-06-07T12:01:29.506698Z","iopub.status.idle":"2023-06-07T12:01:29.517754Z","shell.execute_reply.started":"2023-06-07T12:01:29.506669Z","shell.execute_reply":"2023-06-07T12:01:29.516441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race['program_number'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:29.519315Z","iopub.execute_input":"2023-06-07T12:01:29.519691Z","iopub.status.idle":"2023-06-07T12:01:29.534459Z","shell.execute_reply.started":"2023-06-07T12:01:29.519652Z","shell.execute_reply":"2023-06-07T12:01:29.533091Z"},"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']","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:29.536927Z","iopub.execute_input":"2023-06-07T12:01:29.537332Z","iopub.status.idle":"2023-06-07T12:01:29.560597Z","shell.execute_reply.started":"2023-06-07T12:01:29.53729Z","shell.execute_reply":"2023-06-07T12:01:29.559414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_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":"2023-06-07T12:01:29.562106Z","iopub.execute_input":"2023-06-07T12:01:29.562462Z","iopub.status.idle":"2023-06-07T12:01:29.682936Z","shell.execute_reply.started":"2023-06-07T12:01:29.562431Z","shell.execute_reply":"2023-06-07T12:01:29.681801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":"2023-06-07T12:01:29.685106Z","iopub.execute_input":"2023-06-07T12:01:29.685498Z","iopub.status.idle":"2023-06-07T12:01:38.909944Z","shell.execute_reply.started":"2023-06-07T12:01:29.685465Z","shell.execute_reply":"2023-06-07T12:01:38.908637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_race['sps']","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:38.91209Z","iopub.execute_input":"2023-06-07T12:01:38.912954Z","iopub.status.idle":"2023-06-07T12:01:38.92263Z","shell.execute_reply.started":"2023-06-07T12:01:38.912913Z","shell.execute_reply":"2023-06-07T12:01:38.921691Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:38.924307Z","iopub.execute_input":"2023-06-07T12:01:38.924862Z","iopub.status.idle":"2023-06-07T12:01:38.938787Z","shell.execute_reply.started":"2023-06-07T12:01:38.924828Z","shell.execute_reply":"2023-06-07T12:01:38.937459Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:38.940268Z","iopub.execute_input":"2023-06-07T12:01:38.94061Z","iopub.status.idle":"2023-06-07T12:01:38.953897Z","shell.execute_reply.started":"2023-06-07T12:01:38.940581Z","shell.execute_reply":"2023-06-07T12:01:38.952928Z"},"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":{"execution":{"iopub.status.busy":"2023-06-07T12:01:38.955239Z","iopub.execute_input":"2023-06-07T12:01:38.955971Z","iopub.status.idle":"2023-06-07T12:01:38.969226Z","shell.execute_reply.started":"2023-06-07T12:01:38.955931Z","shell.execute_reply":"2023-06-07T12:01:38.967752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# here I add all the extra information to each individual horse race\n\ndf_corr = add_info(sprints_pp)\ndf_corr = add_d(df_corr)\ndf_corr = df_corr.dropna()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:01:38.970751Z","iopub.execute_input":"2023-06-07T12:01:38.971311Z","iopub.status.idle":"2023-06-07T12:06:59.662791Z","shell.execute_reply.started":"2023-06-07T12:01:38.971278Z","shell.execute_reply":"2023-06-07T12:06:59.661603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# i commented out because i have the file and just read it now\nsprints_draft1 = draft_all(sprints_pp)\nsprints_draft2 =add_draft(sprints_draft1)\ndf_draft = sprints_draft2.copy()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T12:06:59.664815Z","iopub.execute_input":"2023-06-07T12:06:59.665246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# i commented out because i have the file and just read it now\n\nfinal_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_v21',\n       'delta_v32', 'delta_d_race', 'draft_simple', 'draft_v3'\n    ]].copy()\n\n\n\n      \n\ndf_corr2  = df_corr2[(np.abs(stats.zscore(df_corr2)) < 3).all(axis=1)] \n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"lets s start from here and load the finalcorr file","metadata":{}},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (12,12)\n\nsns.pairplot(data=df_corr2, vars = [ 'weight_carried', 'position_at_finish', 'distance',\n       'pa', 'delta_d', 'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts',  'draft_v3'] )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corrs = df_corr2[[  'weight_carried',  'distance',\n       'delta_d', 'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts',  'draft_v3','pa', 'position_at_finish' ]].corr()\nsns.heatmap(df_corrs, annot=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### First observations from the heatmap 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 correlates best with the draft percentage (called here draft v3) but in the opposite direction that we expected: the more a horse was in the draft position, the larger the position at finish. It looks like being in a draft position for large parts of the race **is not good**  for the outcome of the race. \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**\n> - **delta_v10**, **delta_v21**, and **delta_v32** also correlate well with **pa** and inversely correlate with **positon_st_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":"Let's take a closer look at these graphs and correlations by selecting a smaller number of features. ","metadata":{}},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (10,10)\n\nsns.pairplot(data=df_corr2, vars = [ 'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts',\n        'draft_v3','position_at_finish', \n       'pa'] )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corrs = df_corr2[[   'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts' ,'position_at_finish', \n       'pa' ]].corr()\nsns.heatmap(df_corrs, annot=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Indeed the higher the average projected speed the lower position at finish so we can 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 [REF]","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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (12,12)","metadata":{"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')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our data on drafting does not corroborate previously reported data [REF]. The more a horse is in a drafting position during a race, the higher the **finish position**, that is, the lower the performance. This behaviour is seen also when using the simple drafting definition from [REF] (analysis available upon request). ","metadata":{}},{"cell_type":"code","source":"\nplt.rcParams['font.size'] = 18\nplt.rcParams['axes.linewidth'] = 2\nfig = plt.figure(figsize=(6, 6))\nax = fig.add_axes([0.2,0.2,0.85,0.85]) # axis starts at 0.1, 0.1","metadata":{"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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.rcParams['font.size'] = 18\nplt.rcParams['axes.linewidth'] = 2\nfig = plt.figure(figsize=(6,6))\nax = fig.add_axes([0.2,0.1,0.69,0.80]) # axis starts at 0.1, 0.1\nax.xaxis.set_tick_params(which='major', size=8, width=2, direction='in', top='on')\n#ax.xaxis.set_tick_params(which='minor', size=7, width=2, direction='in', top='on')\nax.yaxis.set_tick_params(which='major', size=8, width=2, direction='in', right='on')\n# ax.yaxis.set_tick_params(which='minor', size=7, width=2, direction='in', right='on')\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.title('Position at finish vs. draft for all short races' , fontdict={'fontsize': 18}, loc='center')\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\draft-vs-position-all.jpg\")\nplt.show()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_p = []\nlist_ap = []\nlist_ap_std = []\nfor i, df in df_corr2.groupby('draft_v3'):\n    average = df['position_at_finish'].mean()\n    std_dev = df['position_at_finish'].std()\n    list_p.append(i)\n    list_ap.append(average)\n    list_ap_std.append(std_dev)\n    \nplt.rcParams['font.size'] = 18\nplt.rcParams['axes.linewidth'] = 2\nfig = plt.figure(figsize=(6,6))\nax = fig.add_axes([0.2,0.1,0.70,0.80]) # axis starts at 0.1, 0.1\nax.xaxis.set_tick_params(which='major', size=8, width=2, direction='in', top='on')\n#ax.xaxis.set_tick_params(which='minor', size=7, width=2, direction='in', top='on')\nax.yaxis.set_tick_params(which='major', size=8, width=2, direction='in', right='on')\n# ax.yaxis.set_tick_params(which='minor', size=7, width=2, direction='in', right='on')\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('the draft')\nplt.ylabel('position at finish')\nplt.title('Draft vs outcome for all short races' , fontdict={'fontsize': 18}, loc='center')\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\draft-vs-position-all.jpg\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(list_p)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(list_ap)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr2.columns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.rcParams['font.size'] = 18\nplt.rcParams['axes.linewidth'] = 2\nfig = plt.figure(figsize=(4, 4))\nax = fig.add_axes([0, 0, 1, 1])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ax.xaxis.set_tick_params(which='major', size=2, width=2, direction='in', top='on')\n#ax.xaxis.set_tick_params(which='minor', size=7, width=2, direction='in', top='on')\nax.yaxis.set_tick_params(which='major', size=0.05, width=2, direction='in', right='on')\n# ax.yaxis.set_tick_params(which='minor', size=7, width=2, direction='in', right='on')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"races_d_ft = final_corr.query('course_type == \"D\" & track_condition == \"FT \"')\ndf_corr_4 = races_d_ft[[ 'weight_carried', 'position_at_finish', 'distance',\n       'pa', 'delta_d', 'sigma_v', 'rel_sigma_v',\n       'delta_v_st', 'delta_v_ts', 'delta_v21',\n       'delta_v32', 'delta_d_race', 'draft_simple', 'draft_v3'\n    ]].copy()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.subplots(tight_layout=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_p = []\nlist_ap = []\nlist_ap_std = []\nfor i, df in df_corr_4.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.rcParams['font.size'] = 18\nplt.rcParams['axes.linewidth'] = 2\nfig = plt.figure(figsize=(6,6))   \n\n# https://stackoverflow.com/questions/19576317/matplotlib-savefig-does-not-save-axes\nax = fig.add_axes([0.2,0.1,0.69,0.8]) # axis starts at 0.1, 0.1\nax.xaxis.set_tick_params(which='major', size=8, width=2, direction='in', top='on')\n#ax.xaxis.set_tick_params(which='minor', size=7, width=2, direction='in', top='on')\nax.yaxis.set_tick_params(which='major', size=8, width=2, direction='in', right='on')\n# ax.yaxis.set_tick_params(which='minor', size=7, width=2, direction='in', right='on')\nplt.scatter(x=list_p, y=list_ap ,marker = \"o\",  linewidths = 1)\nplt.errorbar(list_p, list_ap, yerr=list_ap_std, fmt=\"o\")\nplt.xlabel('position at finish')\nplt.ylabel('the draft')\nplt.title('Position at finish vs. draft for \\n course_type \"D\" & track_condition \"FT\"' , fontdict={'fontsize': 18}, loc='center')\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\draft-vs-position-some.jpg\")\nplt.show()","metadata":{"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":"#### 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","metadata":{}},{"cell_type":"code","source":"final_corr_s = final_corr[['track_id','race_date','program_number','course_type','track_condition', 'race_type']]","metadata":{"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_d_race', 'draft_simple',\n       'draft_v3']]","metadata":{"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":{"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)]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_corr = pd.merge(final_corr_n,final_corr_s,left_index=True, right_index=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_corr['furlongs'] = final_corr['distance']/100","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#### Let's model the average projectes speed **pa** on both categorical and numerical features ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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', 'draft_v3' \\\n                 \n                 ]]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X.index.name = 'index'","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = final_corr['pa']\ny.index.name = 'index'","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X.columns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_features = ['purse', 'weight_carried', 'win_prob', 'distance', 'delta_d',\n       'rel_sigma_v', 'delta_v_st', 'delta_v_ts', 'draft_v3', 'course_type_D',\n       'course_type_I', 'course_type_O', 'course_type_T',\n       'track_condition_FM ', 'track_condition_FT ', 'track_condition_GD ',\n       'track_condition_MY ', 'track_condition_SY ', 'track_condition_YL ',\n       'race_type_ALW', 'race_type_AOC', 'race_type_CLM', 'race_type_MCL',\n       'race_type_MSW', 'race_type_SOC', 'race_type_SST', 'race_type_STK',\n       'race_type_STR', 'race_type_WCL', 'race_type_WMC']","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"regressor=xgb.XGBRegressor(eval_metric='rmsle')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/code/carlmcbrideellis/an-introduction-to-xgboost-regression","metadata":{"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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions = regressor.predict(X_test)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean_absolute_error(y_test, predictions)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_test.std()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import eli5\nfrom eli5.sklearn import PermutationImportance\n\n# Make a small change to the code below to use in this problem. \nperm = PermutationImportance(regressor, random_state=1).fit(X_train, y_train)\n\n# Check your answer\n#q_2.check()\n\n# uncomment the following line to visualize your results\neli5.show_weights(perm, feature_names = X_train.columns.tolist())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perm","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pdpbox import pdp, get_dataset, info_plots","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nfeat_name = 'rel_sigma_v'\npdp_dist = pdp.pdp_isolate(model=regressor, dataset=X_test, model_features=base_features, feature=feat_name)\n\npdp.pdp_plot(pdp_dist, feat_name)\nplt.show()\n'''","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nfeat_name = 'draft_v3'\npdp_dist = pdp.pdp_isolate(model=regressor, dataset=X_test, model_features=base_features, feature=feat_name)\n\npdp.pdp_plot(pdp_dist, feat_name)\nplt.show()\n'''","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nfeat_name = 'delta_v_ts'\npdp_dist = pdp.pdp_isolate(model=regressor, dataset=X_test, model_features=base_features, feature=feat_name)\n\npdp.pdp_plot(pdp_dist, feat_name)\nplt.show()\n'''","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(6,6)) ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor, dataset=X_test, model_features=base_features, features=['rel_sigma_v', 'purse'])\n\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['rel_sigma_v', 'purse'], plot_type='contour', plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for rel_sigma_v and purse',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 84,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\n#plt.show()\nplt.savefig(\"D:\\Documents\\Kaggle\\Derby\\paper\\purse_rel_sigma.png\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor, dataset=X_test, model_features=base_features, features=['distance', 'rel_sigma_v'])\n\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['distance', 'rel_sigma_v'], plot_type='contour',plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for rel_sigma_v and distance',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 24,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\n#plt.show()\nplt.savefig(\"D:\\Documents\\Kaggle\\Derby\\paper\\distance_rel_sigma.png\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor, dataset=X_test, model_features=base_features, features=['delta_v_st', 'delta_v_ts'])\n\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['delta_v_st', 'delta_v_ts'], plot_type='contour',plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for delta_v_st and delta_v_ts',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 24,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\n#plt.show()\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\vst_vts.png\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor, dataset=X_test, model_features=base_features, features=['draft_v3', 'rel_sigma_v'])\n\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=[ 'draft_v3', 'rel_sigma_v'], plot_type='contour', plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for rel_sigma_v and drafting',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 24,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\n#plt.show()\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\draft_rel_sigma.png\")","metadata":{"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 projectes speed, in line with our theoretical calculation and the numerical findings in Annex II. In other words, the horse and jockey must try to minimize speed fluctuations as much as possible. \n> - sitting in a low draft 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 course types and track conditions are of utmost importance in the average speed outcome, indicating that the approximation used in previous work [REF] to put them in the model as random variables is not valid. \n> - the race distance expressed meters is a late 7th in importance regarding the average race speed, again in contradiction to the prevous findings [REF] \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 below, the path length, the ** delta_d** feature has very little impact on the average projected speed. ","metadata":{}},{"cell_type":"markdown","source":"### We use XGBoost to model the final position\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. WE chose the races in such a way that the course type and track conditions are the same. ","metadata":{}},{"cell_type":"code","source":"races_d_ft = final_corr.query('course_type == \"D\" & track_condition == \"FT \"')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"races_d_ft.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"races_d_ft.groupby(['track_id','race_date','race_number']).count()","metadata":{"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',  'draft_v3' \\\n                 \n                 ]]\nX5.index.name = 'index'","metadata":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"regressor5=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(regressor5, param_grid, cv=5).fit(X_train5, y_train5)\n\nprint(\"The best hyperparameters are \",search.best_params_)\nregressor5=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')\nregressor5.fit(X_train5, y_train5)\n\npredictions = regressor5.predict(X_test5)\n\nprint('mae is ' , mean_absolute_error(y_test5, predictions))\n\nprint( 'std ground truth is ', y_test5.std())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Make a small change to the code below to use in this problem. \nperm = PermutationImportance(regressor5, random_state=1).fit(X_train5, y_train5)\n\n# Check your answer\n#q_2.check()\n\n# uncomment the following line to visualize your results\neli5.show_weights(perm, feature_names = X_train5.columns.tolist())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor5, dataset=X_test5, model_features=X_test5.columns.to_list(), features=[ 'delta_v_ts', 'win_prob'])\n\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['delta_v_ts', 'win_prob'], plot_type='contour', plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for rel_sigma_v and purse',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 24,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\n#plt.show()\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\delta_v_ts_win_prob2.png\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor5, dataset=X_test5, model_features=X_test5.columns.to_list(), features=['draft_v3', 'delta_v_ts'])\n\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['draft_v3', 'delta_v_ts'], plot_type='contour', plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for draft_v3 and delta_v_ts',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 24,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\n#plt.show()\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\draft_v3_delta_v_ts2.png\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor5, dataset=X_test5, model_features=X_test5.columns.to_list(), features=['delta_v_st', 'draft_v3'])\n\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=[  'delta_v_st', 'draft_v3'], plot_type='contour', plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for delta_v_st and draft_v3',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 24,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\n#plt.show()\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\delta_v_st_draft_v3_2.png\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inter1  =  pdp.pdp_interact(model=regressor5, dataset=X_test5, model_features=X_test5.columns.to_list(), features=['rel_sigma_v', 'delta_v_st'])\npdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=['rel_sigma_v', 'delta_v_st'], plot_type='contour', plot_params = {\n    # plot title and subtitle\n    'title': 'PDP interact for rel_sigma_v and delta_v_st',\n    'subtitle': None,\n    'title_fontsize': 24,\n    'axes_fontsize': 64,\n    # color for contour line\n    'contour_color':  'white',\n    'font_family': 'Arial',\n    # matplotlib color map for interact plot\n    'cmap': 'viridis',\n    # fill alpha for interact plot\n    'inter_fill_alpha': 0.8,\n    # fontsize for interact plot text\n    'inter_fontsize': 18,\n     'fontsize': 18,\n})\nplt.savefig(r\"D:\\\\Documents\\\\Kaggle\\\\Derby\\\\paper\\\\rel_sigma_v_delta_v_st2.png\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stop here","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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":{"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 position at finish, 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** inverse correlation 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 it happened in our previous analysis. ","metadata":{}},{"cell_type":"markdown","source":"### Can we beat the odds? That is, can we 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. <br>\n","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"races_d_ft.groupby(['track_id','race_date','race_number']).count()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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)\nprint('the number of races is ', len(races_list))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"races_list[2]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"random.shuffle(races_list)\nraces_list[2]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_races = races_list[:int(0.8*len(races_list))]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_races = races_list[int(0.8*len(races_list)):]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(train_races) + len(test_races) == len(races_list)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"set(train_races).intersection(set(test_races))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(test_races)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_races[2]","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_df.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train9 = train_df[['weight_carried', 'win_prob', 'delta_v_st']]\nX_train9.index.name = 'index'","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ny_train9 = train_df['position_at_finish']\ny_train9.index.name = 'index'","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_test9 = test_df[['position_at_finish']]\ny_test9.index.name ='index'","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test9 = test_df[['weight_carried', 'win_prob', 'delta_v_st']]\nX_test9.index.name = 'index'","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"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_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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds = pd.Series(predictions, index = y_test9.index, name='prediction')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"positions_preds = pd.merge(y_test9,preds,left_index=True, right_index=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"positions_preds.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test9.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test9 = pd.merge(test_df, positions_preds, left_index=True,right_index=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test9.head(3)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_all.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"accuracy_score(df_all['position_at_finish_x'], df_all['pred_prob'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"accuracy_score(df_all['position_at_finish_x'], df_all['pred_mod'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"0.2738255033557047/0.22550335570469798","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.columns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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',  \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', \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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds = pd.Series(predictions, index = y_test9.index, name='prediction')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"positions_preds = pd.merge(y_test9,preds,left_index=True, right_index=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test9 = pd.merge(test_df, positions_preds, left_index=True,right_index=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"accuracy_score(df_all['position_at_finish_x'], df_all['pred_prob'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"accuracy_score(df_all['position_at_finish_x'], df_all['pred_mod'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\naqu_sprint_d = sprints_pp.query('distance_id < 710 & course_type == \"D\"')\naqu_600_12 = aqu_sprint_d.query('distance_id == 600 & position_at_finish in [1,2] & track_id == \"AQU\"')\naqu_600_5 = aqu_sprint_d.query('distance_id == 600 & position_at_finish > 4 & track_id == \"AQU\"')\ndf_12 = make_df(aqu_600_12)\ndf_5 = make_df(aqu_600_5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Annex I\nThis is just a plot that compares the previously used \"simple draft\" 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":{"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":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nv3s = []\nfor i in indices:\n    speeds  = random.normal(loc=20, scale=2/i, size=200)\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\n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(1/20)*np.random.randn(1,200)+20","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"indices = 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()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}