{"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":"In this work, I have employed a convolutional neural network architecture for the purpose of training a model. The input data for the network consists of a matrix of dimensions (5160 x 6), which comprises all the sensor readings. Each sensor reading is characterized by six distinct features, namely, x, y, z, time, charge, and auxiliary. These input features are fed into the network during the training process. Additionally, I am currently investigating methods for lossless input translation of coordinates. This [literature](https://arxiv.org/abs/2101.11589) suggests several promising approaches in this regard, which are being explored.","metadata":{}},{"cell_type":"code","source":"import os\nimport warnings\nimport math\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport gc\nos.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'\nwarnings.filterwarnings(\"ignore\", category=UserWarning, module='tensorflow')\nimport tensorflow as tf\nfrom tqdm import tqdm\n\ndef 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    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    scalar_prod = sz1*sz2*(ca1*ca2 + sa1*sa2) + (cz1*cz2)\n    scalar_prod =  np.clip(scalar_prod, -1, 1)\n    return np.average(np.abs(np.arccos(scalar_prod)))\n\ndef convert_to_cartesian(azimuth, zenith):\n    x = np.sin(zenith) * np.cos(azimuth)\n    y = np.sin(zenith) * np.sin(azimuth)\n    z = np.cos(zenith)\n    return x, y, z\n\ndef convert_to_polar(x, y, z):\n    r = math.sqrt(x * x + y * y + z * z)\n    c = -1 if y < 0 else 1\n    azimuth = math.acos(x / math.sqrt(x * x + y * y)) * c\n    zenith = math.acos(z / r)\n    return azimuth, zenith\n\ndef adjust_polar(azimuth, zenith):\n    if azimuth < 0:\n        azimuth += math.pi * 2\n    elif zenith < 0:\n        zenith += math.pi\n    return azimuth, zenith","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensors_df = pd.read_csv('./data/sensor_geometry.csv')\nx_min, x_max, y_min, y_max, z_min, z_max = sensors_df.x.min(), sensors_df.x.max(), sensors_df.y.min(), sensors_df.y.max(), sensors_df.z.min(), sensors_df.z.max()\nsensors_df['x'] = (sensors_df['x'] - x_min) / (x_max - x_min)\nsensors_df['y'] = (sensors_df['y'] - y_min) / (y_max - y_min)\nsensors_df['z'] = (sensors_df['z'] - z_min) / (z_max - z_min)\nmeta_df = pd.read_parquet('./data/train_meta.parquet')\nmeta_df[['batch_id', 'event_id', 'first_pulse_index', 'last_pulse_index']] = meta_df[['batch_id', 'event_id', 'first_pulse_index', 'last_pulse_index']].astype(int)\nmeta_df[['azimuth', 'zenith']] = meta_df[['azimuth', 'zenith']].astype(np.float32)\n\ndf = pd.read_parquet('./data/batch_1.parquet')\ndf['auxiliary'] = df['auxiliary'].astype(int)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_size = 5000\n\nmeta_batch_df = meta_df[meta_df.batch_id == 1]\nmeta_batch_df[['target_x', 'target_y', 'target_z']] = meta_batch_df[['azimuth', 'zenith']].apply(lambda x: convert_to_cartesian(x[0], x[1]), axis=1, result_type='expand')\n\nmatrix_main = np.zeros((sample_size, 5160, 6), dtype=np.float32)\ntarget = np.zeros((sample_size, 3))\n\nfor indx in tqdm(range(0, sample_size, 1)):\n    _, event_id, first_pulse_index, last_pulse_index, azimuth, zenith, target_x, target_y, target_z = meta_batch_df.iloc[indx].values\n    df_study = df.iloc[int(first_pulse_index):int(last_pulse_index)].merge(sensors_df, on='sensor_id')\n    minimum_time = df_study['time'].min()\n    period = df_study['time'].max() - minimum_time\n    df_study['time'] = (df_study['time'] - minimum_time) / (period)\n    df_study.set_index('sensor_id', inplace=True)\n    matrix = np.zeros((5160, 6))\n    matrix[df_study.index.values, :] = df_study.values\n\n    matrix_main[indx, :, :] = matrix\n    target[indx, :] = target_x, target_y, target_z\n\ndel df, meta_df, sensors_df, meta_batch_df, df_study\ngc.collect()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_p, validation_p = 0.92, 0.05\n\nmatrix_train = matrix_main[:int(sample_size * train_p), :, :]\ntarget_train = target[:int(sample_size * train_p), :]\n\nmatrix_validation = matrix_main[int(sample_size * train_p):int(sample_size * (train_p + validation_p)), :, :]\ntarget_validation = target[int(sample_size * train_p):int(sample_size * (train_p + validation_p)), :]\n\nmatrix_test = matrix_main[int(sample_size * (train_p + validation_p)):, :, :]\ntarget_test = target[int(sample_size * (train_p + validation_p)):, :]\n\nmatrix_train.shape, target_train.shape, matrix_validation.shape, target_validation.shape, matrix_test.shape, target_test.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.keras.backend.clear_session()\nmodel = tf.keras.models.Sequential([\n    tf.keras.layers.Conv2D(128, (516, 6), activation='relu', input_shape=(5160, 6, 1)),\n    tf.keras.layers.MaxPooling2D((2, 1)),\n    tf.keras.layers.Flatten(),\n    tf.keras.layers.Dense(128, activation='relu'),\n    tf.keras.layers.Dense(64, activation='relu'),\n    tf.keras.layers.Dense(3, activation='linear')\n])\n\n\nmodel.compile(optimizer='adam', loss='mse', metrics=['mae'])\nearly_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=3, restore_best_weights=True)\nmodel.fit(matrix_train, target_train, epochs=100, batch_size=64, validation_data=(matrix_validation, target_validation), callbacks=[early_stopping])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(model.history.history['loss'], label='train')\nplt.plot(model.history.history['val_loss'], label='validation')\nplt.legend()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred = model.predict(matrix_test)\npred_polar = np.array([convert_to_polar(x, y, z) for x, y, z in pred])\ntarget_test_polar = np.array([convert_to_polar(x, y, z) for x, y, z in target_test])\n\naz_true, zen_true = target_test_polar[:, 0], target_test_polar[:, 1]\naz_pred, zen_pred = pred_polar[:, 0], pred_polar[:, 1]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"angular_dist_score(az_true, zen_true, az_pred, zen_pred)","metadata":{},"execution_count":null,"outputs":[]}]}