{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Derby 2022 - Tool Vincenty formulae\n\n","metadata":{}},{"cell_type":"markdown","source":"* To calculate the coordinates of each horse and the distance traveled, I will use [Vincenty's algorithm](https://en.wikipedia.org/wiki/Vincenty%27s_formulae#:~:text=Vincenty's%20formulae%20are%20two%20related,by%20Thaddeus%20Vincenty%20(1975a).) whose margin of error is +- 3 millimeters. \n\n* In calculating the distance traveled by each horse, it is necessary to know the angle at which that distance is traced and the difference with the optimal angle in each section of the route, otherwise the results would be erroneous. For this purpose, each distance from one coordinate to another will be taken as a vector.\n\n# To make your statistics reliable, use Vincenty's formula.\n\n# I leave it here for you to take, it is easy to apply.\n\n# Below is the function to use this formula instead of the Haversine formula which is not suitable for this competition.\n\n# Happy Competition.","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport os\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-10-02T11:26:52.455996Z","iopub.execute_input":"2022-10-02T11:26:52.457035Z","iopub.status.idle":"2022-10-02T11:26:52.461921Z","shell.execute_reply.started":"2022-10-02T11:26:52.456992Z","shell.execute_reply":"2022-10-02T11:26:52.460739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from math import atan\nfrom math import atan2\nfrom math import cos\nfrom math import radians\nfrom math import sin\nfrom math import sqrt\nfrom math import tan\n\n\nclass Vincenty:\n    def __init__(self, c1, c2):        \n        self.meters = self.calc(c1, c2)                   # output distance in meters\n        self.km = self.meters / 1000                      # output distance in kilometers\n        self.mm = self.meters * 1000                      # output distance in millimeters\n        self.miles = self.meters * 0.000621371            # output distance in miles\n        self.n_miles = self.miles * (6080.20 / 5280)      # output distance in nautical miles\n        self.feet = self.miles * 5280                       # output distance in feet\n        self.inches = self.feet * 12                      # output distance in inches\n        self.yards = self.feet / 3                        # output distance in yards\n        \n        \n    def calc(self, c1 ,c2):\n        a = 6378137.0\n        f = 1 / 298.257223563\n        b = (1 - f) * a\n        \n        lat1, lon1, = c1\n        lat2, lon2, = c2\n        \n        U1 = atan((1 - f) * tan(radians(lat1)))\n        U2 = atan((1 - f) * tan(radians(lat2)))\n        \n        L = radians(lon2 - lon1)\n        Lambda = L\n        \n        # COMMON USED VALUES\n        sin_U1 = sin(U1)\n        cos_U1 = cos(U1)\n        sin_U2 = sin(U2)\n        cos_U2 = cos(U2)\n        \n        # VINCENTY CALCULATING DISTANCE IN METERS\n        \n        for iter in range(200):\n            sin_Lambda = sin(Lambda)\n            cos_Lambda = cos(Lambda) \n            sin_Sigma = sqrt((cos_U2 * sin_Lambda) ** 2 + ((cos_U1 * sin_U2) - (sin_U1 * cos_U2 * cos_Lambda)) ** 2 )\n            cos_Sigma = (sin_U1 * sin_U2) + (cos_U1 * cos_U2 * cos_Lambda)\n            Sigma = atan2(sin_Sigma, cos_Sigma)\n            sin_Alpha = (cos_U1 * cos_U2 * sin_Lambda) / sin_Sigma        \n            cos_Sq_Alpha = 1 - (sin_Alpha ** 2)\n        \n            try:\n                cos2_Sigma_M = cos_Sigma - ((2 * sin_U1 * sin_U2) / cos_Sq_Alpha)\n            except ZeroDivisionError:\n                cos2_Sigma_M = 0.0\n        \n            C = (f / 16) * cos_Sq_Alpha * (4 + (f * (4 - 3 * cos_Sq_Alpha)))            \n            LP = Lambda\n            Lambda = L + (1 - C) * f * sin_Alpha * (Sigma + C * sin_Sigma * (cos2_Sigma_M + C * cos_Sigma * (-1 + 2 * (cos2_Sigma_M ** 2))))\n        \n            if (abs(Lambda - LP)) < 1e-12:\n                break\n        else:\n            return None\n            \n        u_Sq = cos_Sq_Alpha * ((a ** 2 - b ** 2) / (b ** 2))\n        #A = 1 + (u_Sq / 16384) * (4096 + u_Sq * (-768 + u_Sq * (320 - 175 * u_Sq)))\n        #B = (u_Sq / 1024) * (256 + u_Sq * (-128 + u_Sq * (74 - 47 * u_Sq)))\n        k1 = (sqrt(1 + u_Sq) - 1) / (sqrt(1 + u_Sq) + 1)\n        A = (1 + 0.25 * k1 **2) / (1 - k1)\n        B = k1 * (1 - (3/8) * k1 ** 2)\n        d_Sigma = B * sin_Sigma * (cos2_Sigma_M + (B / 4) * (cos_Sigma * (-1 + 2 * cos2_Sigma_M ** 2) - (B / 6) * cos2_Sigma_M * (-3 + 4 * sin_Sigma ** 2) * (-3 + 4 * cos2_Sigma_M ** 2)))\n        s = b * A * (Sigma - d_Sigma)\n        Alpha_1 = atan2(cos_U2 * sin_Sigma, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_Sigma)        \n        Alpha_2 = atan2(cos_U1 * sin_Sigma, - sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_Sigma)\n        \n        return s\n        \n        ","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:26:52.464341Z","iopub.execute_input":"2022-10-02T11:26:52.464678Z","iopub.status.idle":"2022-10-02T11:26:52.486693Z","shell.execute_reply.started":"2022-10-02T11:26:52.464649Z","shell.execute_reply":"2022-10-02T11:26:52.485331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# POPULAR HAVERSINE FORMULAE\n\ndef cal_distance(data):\n    \n    R = 6373.0 ## approximate radius of earth in km\n\n    lat1, lon1, lat2, lon2, = data\n  \n\n    dlon = lon2 - lon1\n    dlat = lat2 - lat1\n    a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2\n    c = 2 * atan2(sqrt(a), sqrt(1-a))\n    distance = R * c*1000\n    \n    return radians(distance)","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:26:52.488685Z","iopub.execute_input":"2022-10-02T11:26:52.489033Z","iopub.status.idle":"2022-10-02T11:26:52.506163Z","shell.execute_reply.started":"2022-10-02T11:26:52.489002Z","shell.execute_reply":"2022-10-02T11:26:52.505205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test Vincenty Formulae","metadata":{}},{"cell_type":"code","source":"# TEST VINCENTY FORMULAE \nprint('Vincenty formulae : ',Vincenty([40.670265648245696, -73.82879452057598], [40.6702786191182, -73.8287890775083]).meters)\n\n\n# TEST HAVERSINE FORMULAE \nprint('Haversine formulae : ', cal_distance([40.670265648245696, -73.82879452057598, 40.6702786191182, -73.8287890775083]))    # 5 CENTIMETERS ERROR","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:26:52.507387Z","iopub.execute_input":"2022-10-02T11:26:52.507812Z","iopub.status.idle":"2022-10-02T11:26:52.525602Z","shell.execute_reply.started":"2022-10-02T11:26:52.507775Z","shell.execute_reply":"2022-10-02T11:26:52.524658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_distance(data):              # USING VINCENTY FORMULAE\n    \n    if (data is None)|(data is np.nan):\n        return 0\n\n    lat1, lon1, lat2, lon2, = data[['latitude', 'longitude', 'latitude2', 'longitude2']].values\n    \n    \n    # *--- GROSS DISTANCE  ------------------------------------------------*\n    try:\n        distance = Vincenty([lat1, lon1], [lat2, lon2]).meters     #Distance traveled\n    except:\n        distance = 0\n    \n        \n    return distance\n\ndef calc_angle(data):\n    \n    if (data is None)|(data is np.nan):\n        return 0\n\n    lat1, lon1, lat2, lon2 = data[['latitude', 'longitude', 'latitude2', 'longitude2']].values\n    #lat1, lon1, lat2, lon2 = data[0:]\n    \n    if (lon2-lon1) != 0:\n        angle = atan((lat2 - lat1)/(lon2-lon1))\n    else:\n        angle = atan((lat2 - lat1)/(0.0001))\n       \n    return angle\n\n        \ndef calc_total_distance(data):\n    \n    diff_angle = abs(data['angle_mean'] - data['angle'])    \n    result = data['distance'] * cos(diff_angle)\n    \n    return result","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:26:52.527583Z","iopub.execute_input":"2022-10-02T11:26:52.528338Z","iopub.status.idle":"2022-10-02T11:26:52.539494Z","shell.execute_reply.started":"2022-10-02T11:26:52.528304Z","shell.execute_reply":"2022-10-02T11:26:52.538579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nyra = pd.read_parquet('../input/nyra-2019-completeparquet')\n","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:26:52.540711Z","iopub.execute_input":"2022-10-02T11:26:52.541208Z","iopub.status.idle":"2022-10-02T11:26:54.746578Z","shell.execute_reply.started":"2022-10-02T11:26:52.541178Z","shell.execute_reply":"2022-10-02T11:26:54.745466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"race_type_dict = {'STK': 'Stakes', \n                  'WCL': 'Waiver Claiming', \n                  'WMC': 'Waiver Maiden Claiming',\n                  'SST': 'Starter Stakes',\n                  'SHP': 'Starter Handicap',\n                  'CLM': 'Claiming', \n                  'STR': 'Starter Allowance',\n                  'AOC': 'Allowance Optionl Claimer',\n                  'SOC': 'Starter Optional Claimer',                  \n                  'MCL': 'Maiden Claiming', \n                  'ALW': 'Allowance', \n                  'MSW': 'Maiden Special Weight'\n                 }\n\ntrack_id_dict = {'AQU': 'Aqueduct', \n                 'BEL': 'Belmont', \n                 'SAR': 'Saratoga'\n                }\ncourse_type_dict = {'M': 'Hurdle', \n                    'D': 'Dirt', \n                    'O': 'Outer turf', \n                    'I': 'Inner turf', \n                    'T': 'turf'\n                   }\ntrack_condition_dict = {'YL ': 'Yielding', \n                        'FM ': 'Firm', \n                        'SY ': 'Sloppy', \n                        'GD ': 'Good', \n                        'FT ': 'Fast', \n                        'MY ': 'Muddy', \n                        'SF ': 'Soft'\n                       }\n\nnyra = nyra.replace({'race_type' : race_type_dict})\nnyra = nyra.replace({'track_id' : track_id_dict})\nnyra = nyra.replace({'course_type' : course_type_dict})\nnyra = nyra.replace({'track_condition' : track_condition_dict})\nnyra['race_distance_meters'] = (nyra['distance_id'] / 100) * 201.168\ndata_next = nyra[['race_horse_id','trakus_index','latitude','longitude']].groupby('race_horse_id').shift(-1)\ndata_next.columns = [f'{c}2' for c in data_next.columns]\nnyra = pd.concat([nyra, data_next], axis=1)\n\n","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:26:54.748792Z","iopub.execute_input":"2022-10-02T11:26:54.749127Z","iopub.status.idle":"2022-10-02T11:27:07.810510Z","shell.execute_reply.started":"2022-10-02T11:26:54.749098Z","shell.execute_reply":"2022-10-02T11:27:07.809492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nyra.head(100)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:27:07.811961Z","iopub.execute_input":"2022-10-02T11:27:07.812736Z","iopub.status.idle":"2022-10-02T11:27:07.844905Z","shell.execute_reply.started":"2022-10-02T11:27:07.812690Z","shell.execute_reply":"2022-10-02T11:27:07.843467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"total_dirty = nyra[(nyra['track_id']=='Aqueduct')                     \n                   & (nyra['course_type']=='Dirt') & (nyra['distance_id'] == 600) & (nyra['race_id'] == 2)]\ntotal_dirty['dtraveled'] = total_dirty.apply(lambda x: calc_distance(x), axis=1)\ntotal_dirty['angle'] = total_dirty.apply(lambda x: calc_angle(x), axis=1)\ntotal_dirty","metadata":{"execution":{"iopub.status.busy":"2022-10-02T11:27:07.846584Z","iopub.execute_input":"2022-10-02T11:27:07.847048Z"},"trusted":true},"execution_count":null,"outputs":[]}]}