{"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":"import gc, random\nfrom time import time\nimport numpy as np\nimport pandas as pd\nfrom sklearn.decomposition import PCA\nfrom sklearn.linear_model import LinearRegression\nfrom numba import jit, njit\nimport matplotlib.pyplot as plt\nt0 = time()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-26T22:44:40.095862Z","iopub.execute_input":"2023-01-26T22:44:40.096352Z","iopub.status.idle":"2023-01-26T22:44:40.103215Z","shell.execute_reply.started":"2023-01-26T22:44:40.096319Z","shell.execute_reply":"2023-01-26T22:44:40.102278Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\nBased on the following notebook:\nhttps://www.kaggle.com/code/shlomoron/icecube-eda-pca-baseline-cv-1-28-lb-1-274\n\nInstead of PCA(1) i use eigenvector corresponding to the largest eigenvalue - this is mathematically the same as PCA(1), but more transparent.\nI made some other minor changes.\n","metadata":{}},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"train_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train_meta.parquet')\ntrain_meta = train_meta.loc[train_meta['batch_id']==1].reset_index(drop=True) # only keep batch 1 for now\ntrain = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_1.parquet')\ngeom = pd.read_csv('/kaggle/input/icecube-neutrinos-in-deep-ice/sensor_geometry.csv')\ntrain = train.merge(geom, on='sensor_id', how='left').drop('sensor_id', axis=1)","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:44:40.104870Z","iopub.execute_input":"2023-01-26T22:44:40.105712Z","iopub.status.idle":"2023-01-26T22:45:01.411793Z","shell.execute_reply.started":"2023-01-26T22:44:40.105663Z","shell.execute_reply":"2023-01-26T22:45:01.410231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Scoring function","metadata":{}},{"cell_type":"code","source":"def angular_dist_score(az_true, zen_true, az_pred, zen_pred):    \n    if not (np.all(np.isfinite(az_true)) and\n            np.all(np.isfinite(zen_true)) and\n            np.all(np.isfinite(az_pred)) and\n            np.all(np.isfinite(zen_pred))):\n        raise ValueError(\"All arguments must be finite\")\n    \n    # pre-compute all sine and cosine values\n    sa1 = np.sin(az_true)\n    ca1 = np.cos(az_true)\n    sz1 = np.sin(zen_true)\n    cz1 = np.cos(zen_true)\n    \n    sa2 = np.sin(az_pred)\n    ca2 = np.cos(az_pred)\n    sz2 = np.sin(zen_pred)\n    cz2 = np.cos(zen_pred)\n    \n    # scalar product of the two cartesian vectors (x = sz*ca, y = sz*sa, z = cz)\n    scalar_prod = sz1*sz2*(ca1*ca2 + sa1*sa2) + (cz1*cz2)\n    \n    # scalar product of two unit vectors is always between -1 and 1, this is against nummerical instability\n    # that might otherwise occure from the finite precision of the sine and cosine functions\n    scalar_prod =  np.clip(scalar_prod, -1, 1)\n    \n    # convert back to an angle (in radian)\n    return np.abs(np.arccos(scalar_prod))","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:01.413961Z","iopub.execute_input":"2023-01-26T22:45:01.415031Z","iopub.status.idle":"2023-01-26T22:45:01.425588Z","shell.execute_reply.started":"2023-01-26T22:45:01.414975Z","shell.execute_reply":"2023-01-26T22:45:01.424277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Angle compute logic","metadata":{}},{"cell_type":"code","source":"def compute_angle(xyz, t, charge, aux):\n    # only keep auxiliary = 0\n    idx = (aux == 0)\n    if idx.sum() > 3 and xyz[idx,0].std() + xyz[idx,1].std() > 0.1: # use only aux=0 if enough data and multiple columns\n        xyz     = xyz[idx]\n        t       = t[idx]\n        charge  = charge[idx]\n        aux     = aux[idx]\n        \n    # get largest eigenvector\n    c = np.cov(xyz.T, aweights=((1-aux)*50+1) * charge**2)\n    w, v = np.linalg.eig(c)\n    ii = w.argmax()\n    vector = v[:,ii]\n    \n    # flip if time is anti-correlated with direction\n    d = (xyz * vector.reshape(1, -1)).sum(axis=1)\n    if np.cov(t, d)[0,1] > 0:\n        vector = -vector\n\n    # angles\n    zenith  = np.arccos(vector[2])\n    azimuth = np.arctan2(vector[1], vector[0])\n    if azimuth < 0:\n        azimuth = azimuth + 2 * np.pi\n    w.sort()\n    # scale w so that max=1\n    w = w / w[2]\n    return azimuth, zenith, w[0], w[1], w[2]","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:01.428419Z","iopub.execute_input":"2023-01-26T22:45:01.428870Z","iopub.status.idle":"2023-01-26T22:45:01.445027Z","shell.execute_reply.started":"2023-01-26T22:45:01.428834Z","shell.execute_reply":"2023-01-26T22:45:01.443080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loop over events","metadata":{}},{"cell_type":"code","source":"# loop over events\nN = 10000\nresults = np.zeros([N, 9], dtype=np.float)\nresults[:,0] = train_meta['event_id'].to_numpy()[:N]\nresults[:,1] = train_meta['azimuth'].to_numpy()[:N]\nresults[:,2] = train_meta['zenith'].to_numpy()[:N]\n\nei1 = train_meta['first_pulse_index'].to_numpy()\nei2 = train_meta['last_pulse_index'].to_numpy()\n\nxyz     = train[['x','y','z']].to_numpy()\nt       = train['time'].to_numpy()\ncharge  = train['charge'].to_numpy()\naux     = train['auxiliary'].to_numpy()\ndel train, train_meta\ngc.collect()\nfor i in range(N): # small subset only, for now\n    i1, i2 = ei1[i], ei2[i]\n    results[i,3:8] = compute_angle(xyz[i1:i2,:], t[i1:i2], charge[i1:i2], aux[i1:i2])\nresults[:,8] = angular_dist_score(results[:,1], results[:,2], results[:,3], results[:,4])\nprint('CV=', np.round(results[:,8].mean(), 3), int(time() - t0), 'sec')","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:01.447152Z","iopub.execute_input":"2023-01-26T22:45:01.447809Z","iopub.status.idle":"2023-01-26T22:45:06.423840Z","shell.execute_reply.started":"2023-01-26T22:45:01.447697Z","shell.execute_reply":"2023-01-26T22:45:06.422552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plots","metadata":{}},{"cell_type":"code","source":"df = pd.DataFrame(np.round(results,3))\ndf.columns = ['event_id','a_t','z_t','a_p','z_p','w1','w2','w3','e']\n# error\nh = plt.hist(df['e'], bins=100)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:06.425103Z","iopub.execute_input":"2023-01-26T22:45:06.425430Z","iopub.status.idle":"2023-01-26T22:45:06.817413Z","shell.execute_reply.started":"2023-01-26T22:45:06.425397Z","shell.execute_reply":"2023-01-26T22:45:06.816474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# act vs pred - a\np = plt.plot(df['a_t'],df['a_p'],'.')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:06.818502Z","iopub.execute_input":"2023-01-26T22:45:06.819811Z","iopub.status.idle":"2023-01-26T22:45:07.062951Z","shell.execute_reply.started":"2023-01-26T22:45:06.819770Z","shell.execute_reply":"2023-01-26T22:45:07.061923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# act vs pred - z\np = plt.plot(df['z_t'],df['z_p'],'.')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:07.064308Z","iopub.execute_input":"2023-01-26T22:45:07.065348Z","iopub.status.idle":"2023-01-26T22:45:07.278953Z","shell.execute_reply.started":"2023-01-26T22:45:07.065309Z","shell.execute_reply":"2023-01-26T22:45:07.277653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# azimuth: true and predicted\nh1 = plt.hist(df['a_t'], bins=200, alpha=0.5)\nh2 = plt.hist(df['a_p'], bins=200, alpha=0.5)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:07.281616Z","iopub.execute_input":"2023-01-26T22:45:07.282827Z","iopub.status.idle":"2023-01-26T22:45:08.111586Z","shell.execute_reply.started":"2023-01-26T22:45:07.282771Z","shell.execute_reply":"2023-01-26T22:45:08.110392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"comments on the above histogram of true and predicted azimuth:\n    true is uniformly distributes over 0 to 2 Pi range, as expected\n    predicted has 6 large peaks, and 6 small peaks. Large peaks correspond to events that only trigger 2 vertical columns, so for them azimuth can only take 6 possible values. This is just a property of the data and i don't think we can do anything about it. Impact on error is relatively minor.","metadata":{}},{"cell_type":"code","source":"# zenith: true and predicted\nh1 = plt.hist(df['z_t'], bins=200, alpha=0.5)\nh2 = plt.hist(df['z_p'], bins=200, alpha=0.5)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:08.112889Z","iopub.execute_input":"2023-01-26T22:45:08.113242Z","iopub.status.idle":"2023-01-26T22:45:08.960158Z","shell.execute_reply.started":"2023-01-26T22:45:08.113209Z","shell.execute_reply":"2023-01-26T22:45:08.958618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"comments on the above histogram of true and predicted zenith:\n    true is distributed as sin(x), as expected.\n    predicted: for some reason it is skewed towards lower values - up vs down direction. This appears to have a large impact on error, and so far i cannot figure out what is causing this pattern. If you have an explanation, please put it in comments.","metadata":{}}]}