{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":7878,"databundleVersionId":46689}],"dockerImageVersionId":31328,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Reconstrucción de trayectorias con DBSCAN (Aprendizaje No Supervisado)\nCargando todas las librerías que se emplearán.","metadata":{}},{"cell_type":"code","source":"!pip install git+https://github.com/LAL/trackml-library.git\nimport numpy as np \nimport pandas as pd \nimport zipfile # Descomprimir archivos .zip\nfrom sklearn.cluster import DBSCAN\nfrom sklearn.preprocessing import StandardScaler\nfrom trackml.score import score_event\n# import os  # Solo lo usé para obtener las rutas hazta los archivos .zip\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import TensorDataset, DataLoader\nfrom sklearn.model_selection import train_test_split\nimport time\nimport matplotlib.pyplot as plt\nimport matplotlib.cm as cm\nfrom sklearn.metrics import (confusion_matrix,classification_report,ConfusionMatrixDisplay)\nfrom sklearn.metrics import roc_curve, auc","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T22:44:07.689923Z","iopub.execute_input":"2026-07-08T22:44:07.690400Z","iopub.status.idle":"2026-07-08T22:44:19.975471Z","shell.execute_reply.started":"2026-07-08T22:44:07.690365Z","shell.execute_reply":"2026-07-08T22:44:19.974784Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"'''\n# Obtenemos las rutas a todos los archivos del reto\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n'''\n\n# Nos quedamos con la ruta hacia el paquete de entrenamiento\n# Nos quedaremos con los datos para un solo evento\nruta_train1 = '/kaggle/input/competitions/trackml-particle-identification/train_1.zip'\nruta_detectors = '/kaggle/input/competitions/trackml-particle-identification/detectors.zip'\n\n# Abrimos el archivo .zip en modo lectura\nzipi = zipfile.ZipFile(ruta_train1, 'r')\n\n# Elegimos que evento queremos tratar\nevento = 'event000001001'\n# print(f\"Cargando archivos para el prefijo: {evento}\")\n\n\n# Leemos cada archivo .csv correspondiente a ese evento\nhits = pd.read_csv(zipi.open(f'train_1/{evento}-hits.csv'))\ncells = pd.read_csv(zipi.open(f'train_1/{evento}-cells.csv'))\nparticles = pd.read_csv(zipi.open(f'train_1/{evento}-particles.csv'))\ntruth = pd.read_csv(zipi.open(f'train_1/{evento}-truth.csv'))\n\n# Cerramos el archivo .zip\nzipi.close()\n\n# Abrimos el de detectores y obtenemos sus datos\nzd = zipfile.ZipFile(ruta_detectors, \"r\")\ndetectors = pd.read_csv(zd.open(\"detectors.csv\"))\nzd.close()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T22:44:21.336994Z","iopub.execute_input":"2026-07-08T22:44:21.337816Z","iopub.status.idle":"2026-07-08T22:44:21.978668Z","shell.execute_reply.started":"2026-07-08T22:44:21.337780Z","shell.execute_reply":"2026-07-08T22:44:21.977977Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Definimos una función para ver que partículas son cercanas y etiquetarlas en un grupo.\ndef etiqueta_grupos(hits, eps=0.015):\n    '''\n    Primero extraemos las coordenadas espaciales donde ocurrió el impacto del archivo de hits.\n    Como lo que tenemos es un detector cilindrico, lo mejor es trabajar en coordenadas cilíndricas.\n    Definimos el radio, el ángulo y una magnitud llamada z_rel. Esta es porque como tenemos una \n    especie de hélice cónica, el avance en el eje z entre el radio de la helice es una magnitud constante.\n    '''\n    \n    # Obtenemos las coordenadas y obtenemos el radio y el ángulo correspondiente\n    x = hits['x']\n    y = hits['y']\n    z = hits['z']\n    r = np.sqrt(x**2+y**2)\n    phi = np.arctan2(y,x) \n    theta = np.arctan2(r,z)\n    z_rel = z/(r+10**(-6))\n    eta = -np.log(np.tan(theta/2))\n    # Agrupamos las variables en coordenadas cilindricas en un array y las normalizaremos\n    X = np.column_stack([np.cos(phi), np.sin(phi) ,z_rel, eta])\n    # X = np.column_stack([theta,z_rel])\n    X_scaled = StandardScaler().fit_transform(X)\n\n    # Creamos los clusters, es decir, los grupos de partículas vecinas\n    # Tomamos el algoritmo DBSCAN para saber cuantas partículas hay dentro de la esfera de radio eps\n    # Exigimos que a partir de una particula en el interior de ese radio sea ya un cluster (min_sample = 1)\n    # El algoritmo de optimización usado será el 'kd_tree'\n    clusters = DBSCAN(eps=eps, min_samples=1, algorithm='kd_tree')\n    etiquetas = clusters.fit_predict(X_scaled) # Nos devuelve las etiquetas del grupo al que pertenece cada partícula\n    return etiquetas # Hay etiquetas con -1, significa que son ruido. No pertenecen a ningún grupo. Las eliminaré después.\n \n# Añado a mi hits.cvs una columna extra con estas etiquetas\nhits['etiqueta_grupo'] = etiqueta_grupos(hits)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T22:44:27.394303Z","iopub.execute_input":"2026-07-08T22:44:27.395117Z","iopub.status.idle":"2026-07-08T22:44:28.095805Z","shell.execute_reply.started":"2026-07-08T22:44:27.395037Z","shell.execute_reply":"2026-07-08T22:44:28.095175Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"'''\nPara hacer una comparativa entre las etiquetas de grupo dadas por el DBSCAN y el archivo de truth.csv\nañadimos otra columna al hits.csv que sea el 'particle_id' del archivo truth.csv\n'''\n\ncomparativa = hits[['hit_id', 'etiqueta_grupo']].merge(truth[['hit_id', 'particle_id']], on='hit_id')\n# print(comparativa.head(50))\n\n# Eliminamos el ruido\ncomparativa_sin_ruido = comparativa[comparativa['etiqueta_grupo']!=-1]\n\n# Ahora tengo que saber si las particulas de un mismo grupo son la misma o diferentes\n# Uso .groupby para ver cuantas partículas diferentes hay en cada grupo y .nunique para saber si son únicas\nparticulas_por_grupo = comparativa_sin_ruido.groupby('etiqueta_grupo')['particle_id'].nunique()\n\n# El total de grupos (de trayectorias) será:\ntotal_grupos = len(particulas_por_grupo)\nprint(f\"Total de trayectorias encontradas: {total_grupos}\")\n\n# El número de grupos con una sola partícula en ellos, es decir, un seguimiento perfecto de la trayectoria será:\ntracks_perfectos = (particulas_por_grupo == 1).sum() # Sumamos todas las veces que en un grupo hay una sola partícula\nprint(f\"Trayectorias perfectas de 1 sola partícula: {tracks_perfectos}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T21:52:45.039753Z","iopub.execute_input":"2026-07-08T21:52:45.040520Z","iopub.status.idle":"2026-07-08T21:52:45.070145Z","shell.execute_reply.started":"2026-07-08T21:52:45.040488Z","shell.execute_reply":"2026-07-08T21:52:45.069479Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Creo una función para crear el conjunto que usaremos para realizar obtener el score.\n# Se necesita un Data Frame de 3 columnas, el identificador del evento, el identificador del hit y el identificador de la trayectoria.\ndef fun_submission(hits, labels, event_id=0):\n    sub_data = np.column_stack(([event_id]*len(hits), hits.hit_id.values, labels))\n    submission = pd.DataFrame(data=sub_data, columns=[\"event_id\", \"hit_id\", \"track_id\"]).astype(int)\n    return submission","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T21:24:26.168446Z","iopub.execute_input":"2026-07-08T21:24:26.168999Z","iopub.status.idle":"2026-07-08T21:24:26.175106Z","shell.execute_reply.started":"2026-07-08T21:24:26.168967Z","shell.execute_reply":"2026-07-08T21:24:26.174256Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Representación en función del eps\nepss = np.linspace(0.005, 0.05, 50)\nscores = []\nfor eps in epss:\n    grupos = etiqueta_grupos(hits, eps)\n    submission = fun_submission(hits, grupos)\n    score = score_event(truth, submission[['hit_id', 'track_id']])\n    scores.append(score)\n\nidx_max = np.argmax(scores) # máximo de los scores\neps_opt = epss[idx_max]\nscore_max = scores[idx_max]\n\nprint(f\"Mejor eps   = {eps_opt:.3f}\")\nprint(f\"Mejor score = {score_max*100:.2f}%\")\n\n\n\nplt.figure(figsize=(8,5))\nplt.plot(epss,scores,\"-o\",linewidth=2,markersize=6)\n\nplt.xlabel(r\"$\\varepsilon$\", fontsize=14)\nplt.ylabel(\"Score TrackML\", fontsize=14)\nplt.title(\"Influencia del parámetro $\\\\varepsilon$ de DBSCAN\", fontsize=15)\nplt.grid()\nplt.tight_layout()\nplt.scatter(eps_opt,score_max,color=\"red\",s=90,zorder=3,label=f\"Máximo (ε={eps_opt:.3f}, Score={score_max*100:.2f}%)\")\nplt.legend()\nplt.savefig(\"score_vs_eps.png\", dpi=600)\nplt.savefig(\"score_vs_eps.pdf\")\n\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T21:24:28.306564Z","iopub.execute_input":"2026-07-08T21:24:28.307250Z","iopub.status.idle":"2026-07-08T21:25:15.307218Z","shell.execute_reply.started":"2026-07-08T21:24:28.307217Z","shell.execute_reply":"2026-07-08T21:25:15.306455Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Reconstrucción con aprendizaje supervisado (DNN)\nAqui comienza la segunda parte del codigo que consiste en el entrenamiento de una red densa con la capa de entrada, la de salida y 4 capas ocultas. La red es 37-800-400-400-200-1.\n","metadata":{}},{"cell_type":"code","source":"# DNN\nDBS = True\nvar = 18\nif DBS:\n    dim = var*2 +1\nelse:\n    dim = var*2\nclass tracker(nn.Module):\n    def __init__(self, input_dim=dim):\n        super().__init__()\n\n        self.network = nn.Sequential(\n            nn.Linear(input_dim, 800),\n            nn.SELU(),\n\n            nn.Linear(800, 400),\n            nn.SELU(),\n\n            nn.Linear(400, 400),\n            nn.SELU(),\n\n            nn.Linear(400, 400),\n            nn.SELU(),\n\n            nn.Linear(400, 200),\n            nn.SELU(),\n\n            nn.Linear(200, 1),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x):\n        return self.network(x)\n\n# Ponemos el modelo en la GPU y como tenemos dos, empleando las dos GPU T4 trabajamos en paralelo.\ndevice = torch.device(\"cuda\")\n\nmodelo = tracker().to(device)  # device = cuda:0\n\nif torch.cuda.device_count() > 1:\n    modelo = nn.DataParallel(modelo, device_ids=[0, 1])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T22:44:40.276714Z","iopub.execute_input":"2026-07-08T22:44:40.277324Z","iopub.status.idle":"2026-07-08T22:44:40.909831Z","shell.execute_reply.started":"2026-07-08T22:44:40.277290Z","shell.execute_reply":"2026-07-08T22:44:40.908910Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def features_hits(hits, cells, detectors):\n    grupo = etiqueta_grupos(hits) # Con el DBSCAN asignamos una etiqueta a cada hit\n    hits = hits.merge(detectors, on=[\"volume_id\",\"layer_id\",\"module_id\"], how=\"left\") # Combinamos hits con el archivo de detectors\n    hit_cells = cells.groupby('hit_id').size() # Contamos el número de celdas sobre las que se deposita energía\n    hit_value = cells.groupby('hit_id')['value'].sum() # Obtenemos el total de energía depositada por cada partícula\n\n    # Ordenamos los hits con .reindex.\n    hit_cells = hit_cells.reindex(hits.hit_id).fillna(0).values\n    hit_value = hit_value.reindex(hits.hit_id).fillna(0).values\n\n    # Obtengo las variables espaciales\n    x = hits['x'].values/1000\n    y = hits['y'].values/1000\n    z = hits['z'].values/1000\n    r = np.sqrt(x**2 + y**2)\n    phi = np.arctan2(y,x) \n    theta = np.arctan2(r,z)\n    eta = -np.log(np.tan(theta/2))\n    rho = np.sqrt(x*x+y*y+z*z)\n    zr = z/(r+1e-6)\n    rz = r/(np.abs(z)+1e-6)\n    \n    # Coordenadas geométricas del módulo detector donde se registró el impacto.\n    # Aporta información adicional sobre la geometria del detector.\n    # En metros\n    cx = hits.cx.values/1000\n    cy = hits.cy.values/1000\n    cz = hits.cz.values/1000\n\n    # Vector perpendicular al plano del sensor donde impacta la partícula\n    rot_xw = hits.rot_xw.values\n    rot_yw = hits.rot_yw.values\n    rot_zw = hits.rot_zw.values\n    \n    # Las meto en un mismo array ya normalizadas\n    features = np.column_stack([x, y, z, cx, cy, cz, rot_xw, rot_yw, rot_zw, r, rho, np.cos(phi), np.sin(phi), eta, zr, rz, hit_cells/10, hit_value])\n    return features, grupo\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T22:44:43.290173Z","iopub.execute_input":"2026-07-08T22:44:43.290827Z","iopub.status.idle":"2026-07-08T22:44:43.298817Z","shell.execute_reply.started":"2026-07-08T22:44:43.290794Z","shell.execute_reply":"2026-07-08T22:44:43.297933Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"entreno = []\nfor event_id in range(10,17):\n    # Primero obtengo los datos de los archivos para cada evento. Utilizaré 10 eventos para entrenar la red.\n    evento = f'event0000010{event_id:02d}'\n    zipi = zipfile.ZipFile(ruta_train1, 'r')\n    hits = pd.read_csv(zipi.open(f'train_1/{evento}-hits.csv'))\n    cells = pd.read_csv(zipi.open(f'train_1/{evento}-cells.csv'))\n    particles = pd.read_csv(zipi.open(f'train_1/{evento}-particles.csv'))\n    truth = pd.read_csv(zipi.open(f'train_1/{evento}-truth.csv'))\n    zipi.close()\n    \n    features, grupos = features_hits(hits, cells, detectors) # Variables\n    \n    # Obtenemos los particle_id verdaderos y eliminamos el ruido\n    particle_ids = truth.particle_id.unique()\n    particle_ids = particle_ids[particle_ids != 0]\n\n    # Primero. Creo pares de hits que sí pertenezcan a la misma partícula para entrenar la red neuronal en positivo.\n    pares_positivos = [] # Para guardar los hits producidos por las particulas\n    for p in particle_ids:\n        # Extraigo los hits_id que corresponden a la misma partícula\n        indices = truth.loc[truth.particle_id == p,'hit_id'].values - 1 # Con el -1 reindexamos los hits_id para empezar como en los arrays\n        \n        for i in indices:\n            for j in indices:\n                if i != j: # No queremos los pares (i,i) con los mismos índices\n                    pares_positivos.append([i,j])\n                    \n    pares_positivos = np.array(pares_positivos)\n\n    # Creamos un grupo con etiquetas con DBSCAN\n    same_cluster_p = (grupos[pares_positivos[:,0]] == grupos[pares_positivos[:,1]]).astype(np.float32)\n    \n    # Creamos el conjunto de entrenamiento con los pares de hits que sabemos que pertenecen a las mismas partículas\n    # Obtenemos las features del primer índice del par y el segundo. Además añado un array de 1s indicando que sí pertenecen a la misma partícula\n    entreno1 = np.column_stack([features[pares_positivos[:,0]], features[pares_positivos[:,1]], same_cluster_p, np.ones(len(pares_positivos))])\n\n    # Segundo. Creamos los pares para el entrenamiento en negativo. De manera aleatoria, que no pertenezcan a ninguna trayectoria real\n    # Haremos que el tamaño de estos pares sea 3 veces mayor que el del entrenamiento positivo.\n    n_hits = len(hits)\n    size = len(entreno1)*3\n\n    # Vector con el particle_id de cada hit\n    p = truth.particle_id.values\n\n    # Creo indices aleatorios en un array con valores máximos el número de hits de tamaño 3 veces el del entreno 1\n    i = np.random.randint(0,n_hits,size) \n    j = np.random.randint(0,n_hits,size)\n\n    # Filtramos los pares para el entrenamiento negativo. Necesitamos que los particle_id sean distintos para que no pertenezcan a la misma partícula\n    # y si alguno es 0 tambien es negativo el entrenamiento\n    condicion_negativa = ((p[i] == 0)|(p[i] != p[j]))\n    pares_negativos = np.column_stack([i[condicion_negativa],j[condicion_negativa]])\n\n    # Creamos un grupo con etiquetas con DBSCAN\n    same_cluster_n = (grupos[pares_negativos[:,0]] == grupos[pares_negativos[:,1]]).astype(np.float32)\n    \n    # Obtenemos las features del primer índice del par y el segundo. Añado como antes un array de 0s indicando que NO pertenecen a la misma partícula\n    entreno0 = np.column_stack([features[pares_negativos[:,0]], features[pares_negativos[:,1]], same_cluster_n, np.zeros(len(pares_negativos))])\n    entreno.append(entreno1)\n    entreno.append(entreno0)\n\nentreno = np.vstack(entreno) # Junto verticalmente el array de entreno\n#print(entreno)\nnp.random.shuffle(entreno) # Lo mezclo para que no queden todas las filas agrupadas\n\nX = entreno[:,:-1].astype(np.float32) # Elimino la columna de etiquetas de 1 y 0 para entrenar.\ny = entreno[:,-1].astype(np.float32) # Me quedo únicamente con las etiquetas.\n\n#np.save(\"entreno.npy\", entreno)\n#np.savez_compressed(\"/kaggle/working/entreno_comprimido.npz\", entreno=entreno)\nprint(entreno.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T21:53:04.744492Z","iopub.execute_input":"2026-07-08T21:53:04.745043Z","iopub.status.idle":"2026-07-08T21:54:46.686782Z","shell.execute_reply.started":"2026-07-08T21:53:04.745011Z","shell.execute_reply":"2026-07-08T21:54:46.685925Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Creo las muestras de entrenamiento y validación. reservaré solo un 5% para la validación.\nX_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.05, shuffle=True, random_state=42)\n\n# Convierto todo en tensores, directamente en GPU\nX_train = torch.tensor(X_train, dtype=torch.float32).to(device)\ny_train = torch.tensor(y_train, dtype=torch.float32).reshape(-1,1).to(device)\nX_val = torch.tensor(X_val, dtype=torch.float32).to(device)\ny_val = torch.tensor(y_val, dtype=torch.float32).reshape(-1,1).to(device)\n\n# Ya NO usamos TensorDataset ni DataLoader: todo vive en GPU y barajamos con índices\nperdida = nn.BCELoss()\n\n\nhistorial = []\ndef entrenar(lr, epocas, model=modelo, X=X_train, Y=y_train, loss=perdida, batch_size=8000):\n    device_model = next(modelo.parameters()).device  # Dispositivo\n    optimizador = torch.optim.Adam(model.parameters(), lr=lr) #Uso el optimizador Adam\n    # Número total de ejemplos de entrenamiento y el número de mini-batches necesarios.\n    n = X.shape[0]\n    n_batches = (n + batch_size - 1) // batch_size\n    t = []\n    for epoca in range(epocas):\n        ini = time.time() # Contamos el tiempo que tarda en cada época\n        model.train() # Activamos el entrenamiento\n\n        # Acumuladores en GPU\n        loss_acum = 0.0\n        acc_acum = 0.0\n        \n        perm = torch.randperm(n)  # Shufleamos\n        \n        for i in range(0, n, batch_size):\n            idx = perm[i:i+batch_size] # Seleccionamos las muestras del tamaño del mini batch\n\n            # Transferimos los datos a la GPU\n            x_batch = X[idx].to(device_model)\n            y_batch = Y[idx].to(device_model)\n\n            optimizador.zero_grad() # Acumulo grafientes\n            predic = model(x_batch) # Predicciones\n            lost = loss(predic, y_batch) # Pérdida\n            lost.backward() # Back propagation\n            optimizador.step() # Actualización de los pesos\n\n            with torch.no_grad():\n                    loss_acum += lost.item() #Acumulamos la pérdida\n                    pred_bin = (predic >= 0.5).float() # Tomamos como umbral de clasificación por encima del 0.5\n                    acc_acum += (pred_bin == y_batch).float().mean().item()  # Acumulamos la accuracy\n\n        loss_medio = loss_acum / n_batches # Pérdida media\n        acc_media = acc_acum / n_batches #Accuracy media\n\n        # Ahora evaluamos sobre el conjunto de validación\n        model.eval()\n        with torch.no_grad():\n            pred_val = model(X_val.to(device_model))\n            loss_val = loss(pred_val,y_val.to(device_model)).item() # Pérdida del conjunto de validación\n            pred_bin_val = (pred_val >= 0.5).float()\n            acc_val = (pred_bin_val == y_val.to(device_model)).float().mean().item() # Accuracy del conjunto de validación\n\n        end = time.time()\n        dif = end - ini\n        t.append(dif)\n        historial.append([epoca+1, loss_medio, acc_media, loss_val, acc_val])\n        print(f\"LR={lr:.0e} | Epoch {epoca+1}/{epocas} | Loss={loss_medio:.6f} | \"\n              f\"Acc={acc_media:.4f} | ValLoss={loss_val:.6f} | ValAcc={acc_val:.4f} | Time={dif:.2f} s|\")\n    print(f'Tiempo total = {sum(t):.2f} s')\n    return historial","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T21:54:48.683934Z","iopub.execute_input":"2026-07-08T21:54:48.684984Z","iopub.status.idle":"2026-07-08T21:54:56.646148Z","shell.execute_reply.started":"2026-07-08T21:54:48.684944Z","shell.execute_reply":"2026-07-08T21:54:56.644983Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"entrenar(lr=1e-5, epocas=2)\nentrenar(lr=1e-4, epocas=20)\nentrenar(lr=1e-5, epocas=3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T21:54:58.043728Z","iopub.execute_input":"2026-07-08T21:54:58.044109Z","iopub.status.idle":"2026-07-08T22:09:46.688438Z","shell.execute_reply.started":"2026-07-08T21:54:58.044039Z","shell.execute_reply":"2026-07-08T22:09:46.687541Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Guardamos todo\nnp.save(\"historial.npy\", np.array(historial))\ntorch.save(X_train, \"/kaggle/working/X_train.pt\")\ntorch.save(y_train, \"/kaggle/working/y_train.pt\")\ntorch.save(X_val, \"/kaggle/working/X_val.pt\")\ntorch.save(y_val, \"/kaggle/working/y_val.pt\")\ntorch.save(modelo.module.state_dict(), \"/kaggle/working/modelo_inicial.pt\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:47:16.635152Z","iopub.execute_input":"2026-07-08T19:47:16.636023Z","iopub.status.idle":"2026-07-08T19:47:30.585887Z","shell.execute_reply.started":"2026-07-08T19:47:16.635960Z","shell.execute_reply":"2026-07-08T19:47:30.585101Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"hist = np.array(historial)\n\nloss_train = hist[:,1]\nacc_train  = hist[:,2]\nloss_val   = hist[:,3]\nacc_val    = hist[:,4]\n\n# Numeración continua de épocas\nepocas = np.arange(1, len(hist)+1)\n\n# Cambio de entrenamiento\ncambio_1 = 2\ncambio_2 = 20\n\n# LOSS\nplt.figure(figsize=(8,5))\n\nplt.plot(epocas, loss_train,lw=2,label=\"Entrenamiento\")\n\nplt.plot(epocas, loss_val,lw=2,label=\"Validación\")\n\nplt.axvline(cambio_2+0.5,color=\"black\",linestyle=\"--\",label=\"Cambio de entrenamiento\")\nplt.axvline(cambio_1+0.5,color=\"black\",linestyle=\"--\")\nplt.text(cambio_1-1.95,max(loss_train)*0.95,r\"$\\eta=10^{-5}$\",fontsize=11)\nplt.text(cambio_2-6,max(loss_train)*0.95,r\"$\\eta=10^{-4}$\",fontsize=11)\n\nplt.text(cambio_2+1,max(loss_train)*0.95,r\"$\\eta=10^{-5}$\",fontsize=11)\n\nplt.xlabel(\"Época\")\nplt.ylabel(\"Pérdida\")\nplt.title(\"Evolución de la función de pérdida\")\nplt.grid(alpha=0.3)\nplt.legend()\n\nplt.tight_layout()\nplt.savefig(\"loss_training.png\",dpi=500)\nplt.show()\n\n\n\n# ACCURACY\nplt.figure(figsize=(8,5))\n\nplt.plot(epocas,acc_train*100,lw=2,label=\"Entrenamiento\")\n\nplt.plot(epocas,acc_val*100,lw=2,label=\"Validación\")\n\nplt.axvline(cambio_2+0.5,color=\"black\",linestyle=\"--\",label=\"Cambio de entrenamiento\")\nplt.axvline(cambio_1+0.5,color=\"black\",linestyle=\"--\")\nplt.text(cambio_1-1.95,98,r\"$\\eta=10^{-5}$\",fontsize=11)\nplt.text(cambio_2-6,98,r\"$\\eta=10^{-4}$\",fontsize=11)\nplt.text(cambio_2+1,98,r\"$\\eta=10^{-5}$\",fontsize=11)\n\nplt.xlabel(\"Época\")\nplt.ylabel(\"Accuracy (%)\")\nplt.title(\"Evolución de la exactitud\")\nplt.grid(alpha=0.3)\nplt.legend()\n\nplt.tight_layout()\nplt.savefig(\"accuracy_training.png\",dpi=500)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:59:54.433955Z","iopub.execute_input":"2026-07-08T19:59:54.434542Z","iopub.status.idle":"2026-07-08T19:59:56.389113Z","shell.execute_reply.started":"2026-07-08T19:59:54.434507Z","shell.execute_reply":"2026-07-08T19:59:56.388349Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Hisotgrama de clasificación","metadata":{}},{"cell_type":"code","source":"modelo.eval()\n\nwith torch.no_grad():\n    pred = modelo(X_val.to(device)).cpu().numpy().ravel()\ny_real = y_val.cpu().numpy().ravel()\n\nplt.hist(pred[y_real == 1],bins=100,alpha=0.5,label=\"Positivos\", color = \"forestgreen\")\nplt.hist(pred[y_real == 0],bins=100,alpha=0.5,label=\"Negativos\", color = \"crimson\")\n\nplt.legend()\n\nplt.savefig(\"/kaggle/working/histograma.png\",dpi=300,bbox_inches=\"tight\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:49:02.141454Z","iopub.execute_input":"2026-07-08T19:49:02.142236Z","iopub.status.idle":"2026-07-08T19:49:03.466130Z","shell.execute_reply.started":"2026-07-08T19:49:02.142200Z","shell.execute_reply":"2026-07-08T19:49:03.465465Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# MATRIZ DE CONFUSIÓN\ny_pred = (pred > 0.5).astype(int)\n\ncm = confusion_matrix(y_real, y_pred)\nprint(\"Matriz de confusión:\")\nprint(cm)\n\ndisp = ConfusionMatrixDisplay(\n    confusion_matrix=cm,\n    display_labels=[\"Negativo\", \"Positivo\"]\n)\n\nfig, ax = plt.subplots(figsize=(6,6))\ndisp.plot(ax=ax, colorbar=False)\nax.set_title(\"Matriz de confusión\")\n\nplt.tight_layout()\nplt.savefig(\"matriz_confusion.png\", dpi=600, bbox_inches=\"tight\")\nplt.savefig(\"matriz_confusion.pdf\", bbox_inches=\"tight\")\n\nplt.show()\n\nTN = cm[0,0]\nFP = cm[0,1]\nFN = cm[1,0]\nTP = cm[1,1]\n\nprint(f\"TP = {TP}\")\nprint(f\"TN = {TN}\")\nprint(f\"FP = {FP}\")\nprint(f\"FN = {FN}\")\n\naccuracy = (TP+TN)/(TP+TN+FP+FN)\nprecision = TP/(TP+FP)\nrecall = TP/(TP+FN)\n\nprint(f\"\\nAccuracy  = {accuracy*100:.2f}%\")\nprint(f\"Precision = {precision*100:.2f}%\")\nprint(f\"Recall    = {recall*100:.2f}%\")\n\n\n# Para obtener la curva ROC necesitamos las predicciones del modelo\nmodelo.eval()\n\n# Con esta función se realizarán predicciones por batches para obtener la curva ROC\n\ndef predecir(modelo, X, batch_size=50000):\n    pred = []\n    with torch.no_grad():\n        for i in range(0, len(X), batch_size):\n            batch = X[i:i+batch_size].to(device)\n            salida = modelo(batch)\n            pred.append(salida.cpu().numpy())\n    return np.concatenate(pred).ravel()\n\nprob_train = predecir(modelo, X_train)\nprob_val   = predecir(modelo, X_val)\n\ny_train_np = y_train.cpu().numpy().ravel()\ny_val_np   = y_val.cpu().numpy().ravel()\n\nfpr_train, tpr_train, _ = roc_curve(y_train_np, prob_train)\nfpr_val, tpr_val, _ = roc_curve(y_val_np, prob_val)\n\nauc_train = auc(fpr_train, tpr_train)\nauc_val = auc(fpr_val, tpr_val)\n\nplt.figure(figsize=(7,7))\n\nplt.plot(fpr_train,tpr_train,linewidth=2,label=f\"Entrenamiento (AUC = {auc_train:.4f})\")\n\nplt.plot(fpr_val, tpr_val, linewidth=2, label=f\"Validación (AUC = {auc_val:.4f})\")\n\nplt.plot([0,1],[0,1],'k--',label=\"Clasificador aleatorio\")\n\nplt.xlabel(\"False Positive Rate\")\nplt.ylabel(\"True Positive Rate\")\nplt.title(\"Curva ROC\")\n\nplt.grid(alpha=0.3)\nplt.legend()\n\nplt.tight_layout()\nplt.savefig(\"curva_ROC.png\", dpi=600, bbox_inches=\"tight\")\nplt.savefig(\"curva_ROC.pdf\", bbox_inches=\"tight\")\n\nplt.show()\n\nprint(f\"AUC entrenamiento: {auc_train:.4f}\")\nprint(f\"AUC validación:    {auc_val:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T20:12:25.146681Z","iopub.execute_input":"2026-07-08T20:12:25.147718Z","iopub.status.idle":"2026-07-08T20:13:19.362624Z","shell.execute_reply.started":"2026-07-08T20:12:25.147683Z","shell.execute_reply":"2026-07-08T20:13:19.361493Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"modelo.load_state_dict(torch.load(\"/kaggle/input/datasets/manuelsnchezrodrguez/con-dbscan/modelo_inicial.pt\"))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T20:24:15.302259Z","iopub.execute_input":"2026-07-08T20:24:15.308451Z","iopub.status.idle":"2026-07-08T20:24:15.434088Z","shell.execute_reply.started":"2026-07-08T20:24:15.308398Z","shell.execute_reply":"2026-07-08T20:24:15.433208Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Entrenamiento con negativos difíciles","metadata":{}},{"cell_type":"code","source":"Train_hard = []\nmodelo.eval()\n\n# Lo hacemos para 15 eventos\nfor event_id in range(10,25):\n    # print(f\"Evento {event_id}\")\n    evento = f'event0000010{event_id:02d}'\n    zipi = zipfile.ZipFile(ruta_train1,'r')\n    hits = pd.read_csv(zipi.open(f'train_1/{evento}-hits.csv'))\n    cells = pd.read_csv(zipi.open(f'train_1/{evento}-cells.csv'))\n    truth = pd.read_csv(zipi.open(f'train_1/{evento}-truth.csv'))\n    zipi.close()\n\n    features, grupos = features_hits(hits,cells,detectors)\n    \n    # PARES NEGATIVOS ALEATORIOS. Proceso similar al del primer entrenamiento.\n    size=len(entreno1)*3\n    n_hits = len(hits)\n    i = np.random.randint(0,n_hits,size)\n    j = np.random.randint(0,n_hits,size)\n    p = truth.particle_id.values\n    mask = ((p[i] == 0)|(p[i] != p[j]))\n    pares_negativos = np.column_stack([i[mask],j[mask]])\n    same_cluster = (grupos[pares_negativos[:,0]] ==grupos[pares_negativos[:,1]]).astype(np.float32)\n    Train0 = np.column_stack([features[pares_negativos[:,0]], features[pares_negativos[:,1]], same_cluster, np.zeros(len(pares_negativos))])\n    \n    # Predicciones por batches\n    predicciones = []\n    batch_size = 30000\n    with torch.no_grad():\n        for inicio in range(0,len(Train0),batch_size):\n            fin = min(inicio + batch_size,len(Train0))\n            # Creamos el tensor en la CPU de forma normal\n            batch = torch.tensor(Train0[inicio:fin, :-1], dtype=torch.float32)\n            pred_batch = modelo(batch)   # Introducimos el batch en el modelo\n            predicciones.append(pred_batch.cpu().numpy())  # AÑadimos esas predicciones\n            del batch\n            del pred_batch\n    pred = np.concatenate(predicciones).ravel()\n\n    # Los negativos difíciles serán aquellos que se encientren en un rango entre 0.5 y 0.95 de probabilidad.\n    idx_hard = np.where((pred > 0.5) & (pred < 0.95))[0]\n    Train_hard.append(Train0[idx_hard])\n\nTrain_hard = np.vstack(Train_hard)\nprint(\"Entreno original:\", len(entreno))\nprint(\"Negativos difíciles:\", len(Train_hard))\nprint(f\"Fracción: {len(Train_hard)/len(entreno)*100:.2f}%\")\nprint(\"Negativos difíciles:\",Train_hard.shape)\n\n# Los unimos, añadiendo varias veces los train_hard\nentreno_hard_arr = np.vstack([entreno, Train_hard, Train_hard, Train_hard])\nnp.random.shuffle(entreno_hard_arr)\n\nX_hard = torch.tensor(entreno_hard_arr[:,:-1], dtype=torch.float32)\ny_hard = torch.tensor(entreno_hard_arr[:,-1], dtype=torch.float32).reshape(-1,1)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T22:43:44.366721Z","iopub.execute_input":"2026-07-08T22:43:44.367505Z","execution_failed":"2026-07-08T22:43:50.791Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Entrenamiento con negativos difíciles\nentrenar(lr=1e-4,epocas=30,model=modelo, X=X_hard, Y=y_hard, loss=perdida)\nentrenar(lr=1e-5,epocas=10, model=modelo, X=X_hard, Y=y_hard, loss=perdida)\nentrenar(lr=1e-6,epocas=2, model=modelo, X=X_hard, Y=y_hard, loss=perdida)\n\nwith torch.no_grad():\n    pred = modelo(X_val).cpu().numpy().ravel()  # modelo_hard, no modelo\ny_real = y_val.cpu().numpy().ravel()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T22:13:36.184261Z","iopub.execute_input":"2026-07-08T22:13:36.185202Z","iopub.status.idle":"2026-07-08T22:42:47.356325Z","shell.execute_reply.started":"2026-07-08T22:13:36.185161Z","shell.execute_reply":"2026-07-08T22:42:47.355481Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"np.save(\"hard_negatives.npy\", Train_hard)\ntorch.save(modelo.module.state_dict(), \"/kaggle/working/modelo_ha.pt\")\ntorch.save(X_hard, \"/kaggle/working/X_hard.pt\")\ntorch.save(y_hard, \"/kaggle/working/y_hard.pt\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"modelo.load_state_dict(torch.load(\"/kaggle/input/datasets/manuelsnchezrodrguez/con-dbscan/modelo_hard.pt\"))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# HISTOGRAMA DE PREDICCIONES\n\nplt.figure(figsize=(8,5))\n\nplt.hist(pred[y_real == 1],bins=100,alpha=0.5,label=\"Positivos\", color = \"forestgreen\")\nplt.hist(pred[y_real == 0],bins=100,alpha=0.5,label=\"Negativos\", color = \"crimson\")\n\nplt.xlabel(\"Salida de la red\")\nplt.ylabel(\"Número de ejemplos\")\nplt.legend()\nplt.grid(alpha=0.3)\nplt.savefig(\"/kaggle/working/histograma_hard.png\",dpi=300,bbox_inches=\"tight\")\nplt.show()\n\n# MATRIZ DE CONFUSIÓN\n\ny_pred = (pred > 0.5).astype(int)\n\ncm = confusion_matrix(y_real,y_pred)\n\nprint(\"\\nMatriz de confusión:\")\nprint(cm)\n\nprint(\"\\nClassification report:\\n\")\n\nprint(classification_report(y_real,y_pred,digits=4))\n\ndisp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=[\"Negativo\",\"Positivo\"])\ndisp.plot()\nplt.title(\"Matriz de confusión\")\nplt.show()\n\nTN = cm[0,0]\nFP = cm[0,1]\nFN = cm[1,0]\nTP = cm[1,1]\n\nprint(f\"TP = {TP}\")\nprint(f\"TN = {TN}\")\nprint(f\"FP = {FP}\")\nprint(f\"FN = {FN}\")\n\naccuracy = (TP + TN) / (TP + TN + FP + FN)\nprecision = TP / (TP + FP)\nrecall = TP / (TP + FN)\n\nprint(f\"\\nAccuracy  = {accuracy*100:.2f}%\")\nprint(f\"Precision = {precision*100:.2f}%\")\nprint(f\"Recall    = {recall*100:.2f}%\")\n\n# Para obtener la curva ROC necesitamos las predicciones del modelo\nmodelo.eval()\ndef predecir(modelo, X, batch_size=50000):\n    pred = []\n    with torch.no_grad():\n        for i in range(0, len(X), batch_size):\n            batch = X[i:i+batch_size].to(device)\n            salida = modelo(batch)\n            pred.append(salida.cpu().numpy())\n    return np.concatenate(pred).ravel()\nprob_train = predecir(modelo, X_hard)\nprob_val   = predecir(modelo, X_val)\n\ny_train_np = y_hard.cpu().numpy().ravel()\ny_val_np   = y_val.cpu().numpy().ravel()\n\nfpr_train, tpr_train, _ = roc_curve(y_train_np, prob_train)\nfpr_val, tpr_val, _ = roc_curve(y_val_np, prob_val)\n\nauc_train = auc(fpr_train, tpr_train)\nauc_val = auc(fpr_val, tpr_val)\n\nplt.figure(figsize=(7,7))\n\nplt.plot(fpr_train,tpr_train,linewidth=2,label=f\"Entrenamiento (AUC = {auc_train:.4f})\")\n\nplt.plot(fpr_val,tpr_val,linewidth=2,label=f\"Validación (AUC = {auc_val:.4f})\")\n\nplt.plot([0,1],[0,1],'k--',label=\"Clasificador aleatorio\")\n\nplt.xlabel(\"False Positive Rate\")\nplt.ylabel(\"True Positive Rate\")\nplt.title(\"Curva ROC\")\n\nplt.grid(alpha=0.3)\nplt.legend()\n\nplt.tight_layout()\nplt.savefig(\"curva_ROC.png\", dpi=600, bbox_inches=\"tight\")\nplt.savefig(\"curva_ROC.pdf\", bbox_inches=\"tight\")\n\nplt.show()\n\nprint(f\"AUC entrenamiento: {auc_train:.4f}\")\nprint(f\"AUC validación:    {auc_val:.4f}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T21:23:14.013740Z","iopub.execute_input":"2026-07-08T21:23:14.014516Z","execution_failed":"2026-07-08T21:23:44.599Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"modelo.load_state_dict(torch.load(\"/kaggle/input/datasets/manuelsnchezrodrguez/con-dbscan/modelo_hard.pt\"))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T23:17:48.408327Z","iopub.execute_input":"2026-07-08T23:17:48.408791Z","iopub.status.idle":"2026-07-08T23:17:48.505678Z","shell.execute_reply.started":"2026-07-08T23:17:48.408759Z","shell.execute_reply":"2026-07-08T23:17:48.505070Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Evaluamos el evento 1001","metadata":{}},{"cell_type":"code","source":"# Obtendremos las trayectorias para el elveno 1001\nevento = 'event000001001'\nzipi = zipfile.ZipFile(ruta_train1,'r')\nhits = pd.read_csv(zipi.open(f'train_1/{evento}-hits.csv'))\ncells = pd.read_csv(zipi.open(f'train_1/{evento}-cells.csv'))\ntruth = pd.read_csv(zipi.open(f'train_1/{evento}-truth.csv'))\nzipi.close()\n\nfeatures, grupos = features_hits(hits, cells, detectors) # Variables\n# np.save(\"features_event001.npy\", features)\n\n# Contamos los impactos registrados en cada módulo del detector\ncount = hits.groupby(['volume_id','layer_id','module_id'])['hit_id'].count().values \n\n# Obtenemos en qué modulo del detector cae cada impacto\nmodule_id = np.zeros(len(hits),dtype=np.int32)\nfor i in range(len(count)):\n    si = np.sum(count[:i])\n    module_id[si:si+count[i]] = i\n    \nden = 37 # neuronas en la capa de entrada\nen = 18 # número de features\n\n\n\n# Con esta función se calcula, para un hit, la probabilidad de que el resto pertenezcan a al trayectoria\ndef get_predict(hit, thr=0.5):\n    Tx = np.zeros((len(features), den), dtype=np.float32) # Matriz de entrada\n    Tx[:, :en] = np.tile(features[hit], (len(features), 1)) # Añadimos las features del hit\n    Tx[:, en:2*en] = features # La segunda parte son el resto de features del resto de hits\n    Tx[:, 2*en] = (grupos[hit] == grupos).astype(np.float32) # Variable que nos dice que impactos pertenecen al mismo grupo generado con DBSCAN\n\n    # Para cada batch devolvemos la predicción asociada\n    with torch.no_grad():\n        pred = []\n        batch_size = 20000\n        for i in range(0,len(Tx),batch_size):\n            batch = torch.tensor(Tx[i:i+batch_size],dtype=torch.float32).to(device)\n            pred.append(modelo(batch).cpu().numpy())\n\n        pred = np.concatenate(pred).ravel()\n    # Realizamos un modelo de inferencia para las probabilidades llamado Test Time Argumentation TTA\n    # Los índices para los cuales la probabilidad supere el umbral establecido realizaremos un intercambio\n    idx = np.where(pred > thr)[0] \n    if len(idx) > 0:\n        Tx2 = np.zeros((len(idx),den),dtype=np.float32)\n        # Intercambiamos los dos hits\n        Tx2[:, :en] = Tx[idx, en:2*en]\n        Tx2[:, en:2*en] = Tx[idx, :en]\n        \n        # same_cluster permanece igual\n        Tx2[:, 2*en] = Tx[idx, 2*en]\n\n        # Obtenemos las predicciones para los datos ahora intercambiados\n        with torch.no_grad():\n            pred1 = []\n            for i in range(0,len(Tx2),batch_size):\n                batch = torch.tensor(Tx2[i:i+batch_size],dtype=torch.float32).to(device)\n                pred1.append(modelo(batch).cpu().numpy())\n            pred1 = np.concatenate(pred1).ravel()\n            \n        # La probabilidad final sera la media entre las dos\n        pred[idx] = (pred[idx]+ pred1) / 2\n\n    return pred\n\n# FUNCIÓN PARA FORMAR LA TRAYECTORIA\ndef get_path(hit, mask, thr):\n    path = [hit] # Inicio de la trayectoria, solo un hit\n    a = 0\n    \n    while True:\n        c = get_predict(path[-1],thr/2) # Calculamos las probabilidades con respecto al ultimo camino añadido\n        mask = (c > thr) * mask # Eliminamos los candidatos que no cumplen el umbral\n        mask[path[-1]] = 0 # Tambien se elimina ese impacto\n\n        # Ahora donde la probabilidad SÍ supere el umbral\n        cand = np.where(c > thr)[0]\n        if len(cand) > 0:\n            # Impedimos que una misma trayectoria contenga dos impactos registrados en el mismo módulo\n            # (una partícula no puede producir dos impactos a la vez)\n            mask[cand[np.isin(module_id[cand],module_id[path])]] = 0\n        \n        a = (c + a) * mask # Acumulamos la probabilidad\n        \n        # Si ningun candidato supera el criterio salimos del bucle\n        if a.max() < thr*len(path):\n            break\n        # Sino añadimos el candidato al camino\n        path.append(int(a.argmax()))\n\n    return sorted(path)\n\n# Devuelvo los 10 primeros hits\nfor hit in range(10):\n    path = get_path(hit,np.ones(len(truth)),0.9)\n    gt = np.where(truth.particle_id.values==truth.particle_id.iloc[hit])[0]\n\n    print()\n    print(\"hit_id =\", hit+1)\n    print(\"reconstruct :\", path)\n    print(\"ground truth:\", gt.tolist())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-09T00:04:19.426904Z","iopub.execute_input":"2026-07-09T00:04:19.427794Z","iopub.status.idle":"2026-07-09T00:04:30.884667Z","shell.execute_reply.started":"2026-07-09T00:04:19.427759Z","shell.execute_reply":"2026-07-09T00:04:30.883704Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Comienza la reconstrucción ","metadata":{}},{"cell_type":"code","source":"'''\nTomo un valor fijo de hits que voy a usar. Solo trabajaré con los primeros 50000 hits de los 93680 que hay.\nEsto es por un tema de tiempo, ya que para obtener un resultado sobre los 93680 hits puede durar varias horas\nUsando los 50000 primeros solo tarda una hora (que ya es bastante)\n'''\nN = 20000\n\nhits = hits.iloc[:N].reset_index(drop=True)\n\ntruth = truth.iloc[:N].reset_index(drop=True)\n\nhit_ids_validos = hits.hit_id.values\n\ncells = cells[cells.hit_id.isin(hit_ids_validos)].reset_index(drop=True)\n\nprint(\"Hits :\", len(hits))\nprint(\"Truth:\", len(truth))\nprint(\"Cells:\", len(cells))\n\nfeatures, grupos = features_hits(hits,cells, detectors) # Calculo las variables \n\npreds = [] # Creo una lista vacía con las predicciones que añadiremos\n\n# Crearemos dos matrices vacías en las que añadiremos la información de los hits\nTestX1 = np.zeros((len(features), den), dtype=np.float32)\nTestX1[:, :en] = features # Guardamos en la primera mitad la información del primer impacto\nTestX = np.zeros((len(features), den), dtype=np.float32)\nTestX[:, en:2*en] = features # Guardamos en la segunda mitad la información del segundo \n\n\nbatch_size = 4000 # Tamaño de lote empleado en la inferencia\n\n# Recorremos todos los hits\nfor i in range(len(features)-1):\n    if i % 1000 == 0:\n        print(i,\"/\",len(features))\n        \n    # Cada hit sera el primer elemento del doblete\n    TestX[i+1:,:en] = np.tile(features[i],(len(features)-i-1,1))\n    \n    # Se añade además la variable del DBSCAN\n    TestX[i+1:, 2*en] = (grupos[i] == grupos[i+1:]).astype(np.float32)\n    \n    pred_total = []\n    # Obtenemos las predicciones correspondientes. Desactivamos los gradientes porque no entrenamos\n    with torch.no_grad():\n        for k in range(0,len(TestX)-i-1,batch_size):\n            batch = torch.tensor(TestX[i+1+k:min(i+1+k+batch_size,len(TestX))],dtype=torch.float32).to(device) # A la GPU\n            pred_total.append(modelo(batch).cpu().numpy()) # Predicciones\n    pred = np.concatenate(pred_total).ravel() # Unimos todas\n\n    idx = np.where(pred > 0.5)[0] # Solo nos quedamos con aquellas con probabilidad mayor a 0.5\n\n    # Proceso de TTA\n    with torch.no_grad():\n        Tx2 = np.zeros((len(idx), den), dtype=np.float32)\n        \n        # Intercambiamos los dos hits de posición\n        Tx2[:, :en] = TestX[idx+i+1, en:2*en]\n        Tx2[:, en:2*en] = TestX[idx+i+1, :en]\n        # same_cluster permanece igual\n        Tx2[:, 2*en] = TestX[idx+i+1, 2*en]\n        Xtta = torch.tensor(Tx2, dtype=torch.float32).to(device)\n        pred1 = modelo(Xtta).cpu().numpy().ravel()\n\n    pred[idx] = (pred[idx]+pred1) / 2  # Se realiza una media entre las dos predicciones.\n        \n    idx = np.where(pred > 0.5)[0] # Solo nos quedamos con aquellas con una probabilidad de más del 0.5\n    preds.append([idx+i+1,pred[idx]])\n    \npreds.append([np.array([],dtype=np.int64),np.array([],dtype=np.float32)])\n\n# Ahora recorremos el bucle hacia atrás para ver que conexiones hay entre los hits. Para conectarlos\nfor i in range(len(preds)):\n    ii = len(preds)-i-1 # Empezamos por el final\n    \n    # Recorremos los hits que esten conectados con el hit ii\n    for j in range(len(preds[ii][0])):\n        jj = preds[ii][0][j] \n        preds[jj][0] = np.insert(preds[jj][0],0,ii) # Aqui añadimos ese vecino a la conexión que había\n        preds[jj][1] = np.insert(preds[jj][1],0,preds[ii][1][j]) # Aquí añadimos la probabilidad asociada a esa conexion\n\n# Los guardamos\npreds = np.array(preds,dtype=object)\nprint(\"preds construido\")\nnp.save(\"preds_event001.npy\", preds)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-09T00:04:37.951664Z","iopub.execute_input":"2026-07-09T00:04:37.952743Z","iopub.status.idle":"2026-07-09T00:15:44.311203Z","shell.execute_reply.started":"2026-07-09T00:04:37.952694Z","shell.execute_reply":"2026-07-09T00:15:44.310164Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"preds = []  # Lista vacía donde se guardarán las predicciones\nTestX = np.zeros((len(features), den), dtype=np.float32) # Matriz vacía para añadir las entradas de la red\nTestX[:, en:2*en] = features # En la segunda mitad se añaden las variables del segundo hit\n\nbatch_size = 50000 # Tamaño del batch\n\nfor i in range(len(features)-1):\n    if i % 1000 == 0:\n        print(i, \"/\", len(features))\n        \n    # Cada hit sera el primer elemento del doblete\n    TestX[i+1:,:en] = np.tile(features[i],(len(TestX)-i-1,1))\n    # Rellenamos los candidatos con los hits posteriores\n    TestX[i+1:, en:2*en] = features[i+1:]\n    \n    # Variable del DBSCAN\n    TestX[i+1:, 2*en] = (grupos[i] == grupos[i+1:]).astype(np.float32)\n\n    # Predicciones, como antes\n    pred_total = []\n    with torch.no_grad():\n        for k in range(0,len(TestX)-i-1,batch_size):\n            batch = torch.tensor(TestX[i+1+k :min(i+1+k+batch_size,len(TestX))],dtype=torch.float32).to(device)\n            pred_total.append(modelo(batch).cpu().numpy())\n    pred = np.concatenate(pred_total).ravel()\n\n    # Primero observamos cuales de ellos tienen una probabilidad superior al 0.2\n    idx = np.where(pred > 0.2)[0]\n\n    # Si los hay obtenemos probabilidades y sino, no.\n    if len(idx) > 0:\n        Tx2 = np.zeros((len(idx), den), dtype=np.float32)\n\n        # Intercambiamos hit1 y hit2\n        Tx2[:, :en] = TestX[idx+i+1, en:2*en]\n        Tx2[:, en:2*en] = TestX[idx+i+1, :en]\n        \n        # same_cluster permanece igual\n        Tx2[:, 2*en] = TestX[idx+i+1, 2*en]\n\n        # Predicciones\n        pred1_total = []\n        with torch.no_grad():\n            for k in range(0, len(idx), batch_size):\n                batch = torch.tensor(Tx2[k:k+batch_size],dtype=torch.float32).to(device)\n                pred1_total.append(modelo(batch).cpu().numpy())\n        \n        pred1 = np.concatenate(pred1_total).ravel()\n        \n    # Ahora comprobamos aquellas con probabilidad mayor a un 0.5\n    idx = np.where(pred > 0.5)[0]\n    preds.append([idx+i+1,pred[idx]])\npreds.append([np.array([],dtype=np.int64),np.array([],dtype=np.float32)])\n\n# Repetimos el proceso de conectar los hits recorriendo el vector de preds en sentido contrario\nfor i in range(len(preds)):\n    ii = len(preds)-i-1\n    for j in range(len(preds[ii][0])):\n        jj = preds[ii][0][j]\n        preds[jj][0] = np.insert(preds[jj][0],0,ii)\n        preds[jj][1] = np.insert(preds[jj][1],0,preds[ii][1][j])\n\npreds = np.array(preds,dtype=object)\nprint()\nprint(\"preds construido\")\n\n# La función de predicciones 2 nos devuelve un vector con las probabilidades asociadas a cada hit que se introduzca\ndef get_predict2(p):\n    c = np.zeros(len(preds),dtype=np.float32)\n    c[preds[p][0]] = preds[p][1]\n    return c\n\n\ndef get_path2(hit, mask, thr):\n    path = [hit] # Inicio de la trayectoria, solo un hit\n    a = 0\n    while True:\n        c = get_predict2(path[-1]) # Probabilidades del ultimo hit añadido a la trayectoria\n        mask = (c > thr) * mask # Solo los candidatos que superen el umbral pasan\n        mask[path[-1]] = 0 # No se vulve a seleccionar ese hit\n        cand = np.where(c > thr)[0] # Los candidatos serán aquellos que superen el umbral\n\n        # ELiminamos los hits pertenecientes al mismo módulo\n        if len(cand) > 0:\n            mask[cand[np.isin(module_id[cand],module_id[path])]] = 0\n        \n        a = (c + a) * mask # Acumulamos las probabilidades\n        # Si ningun candidato alcanza la puntuación mínima exigida se para\n        if a.max() < thr * len(path):\n            break\n            \n        path.append(int(a.argmax())) # se añaden al camino aquellos que tengan la mayor probabilidad\n    return path\n\n\ntracks_all = [] # Lista para almacenar las trayectorias\n\nthr = 0.85 # Umbral\nmulti = True # Para activar el método multipaso\n\n# Creamos un bucle con cada iteracion para un hit\n# MULTIPASO\nfor hit in range(len(preds)):\n    if hit % 1000 == 0:\n        print(hit, \"/\", len(preds))\n\n    m = np.ones(len(truth))\n    path = get_path2(hit,m,thr) # Pasamos el hit por la función que reconstruye trayectorias\n\n    # Sí la longitud del camino es mayor a 1 y el metodo multipaso esta activo se elimina el segundo hit para ver si hay otra mas larga\n    if multi and len(path) > 1:\n        m[path[1]] = 0\n        path2 = get_path2(hit,m,thr)\n\n        # Si la longitud del primero menor entonces los igualamos y volvemos a eliminar el segundo se repite el proceso\n        # El camino más largo suele ser el mejor\n        if len(path) < len(path2):\n            path = path2\n            m[path[1]] = 0\n            path2 = get_path2(hit,m,thr) \n            # Probamos otro camino diferente. Si ahora el siguiente es más largo nos quedamos con ese\n            if len(path) < len(path2):\n                path = path2\n                \n        # Si la longitud del segundo camino es mayor a 1 eliminamos el segundo hit de la segunda trayectoria\n        # Mantenemos el de la primera y volvemos a calcular. SI ahora obtenemos una mejor trayectoria nos quedamos con ese\n        elif len(path2) > 1:\n            m[path[1]] = 1\n            m[path2[1]] = 0\n            path2 = get_path2(hit,m,thr)\n            if len(path) < len(path2):\n                path = path2\n\n    tracks_all.append(path) # Añadimos los caminos a la lista de trayectorias\nprint(\"tracks_all construido\")\ntracks = np.zeros(len(hits), dtype=np.int32)\n\nfor track_id, path in enumerate(tracks_all, start=1):\n    for hit in path:\n        tracks[hit] = track_id # Damos a cada hit en la trayectoria su propio identificador de trayectoria\ntracks = np.array(tracks, dtype=np.int32)\nnp.save(\"tracks_event001.npy\",np.array(tracks, dtype=object)) # Guardamos\nsubmission = pd.DataFrame({\"hit_id\":truth.hit_id.values,\"track_id\":tracks.astype(int)})\n\nscore = score_event(truth,submission)\nprint(\"-------------------------------\")\nprint(f\"Score ={score*100:.2f}%\")\n# ----------------------------------------------------------\n# ESTADÍSTICAS\n# ----------------------------------------------------------\n\nn_conexiones = sum(len(p[0]) for p in preds)\n\nlongitudes = np.array([len(p[0]) for p in preds])\n\nprint(\"Conexiones totales:\",n_conexiones)\nprint(\"Media conexiones:\",n_conexiones/len(preds))\nprint(\"Máximo conexiones:\",longitudes.max())\n\n# Guardamos la reconstrucción buena\ntracks_mp = tracks.copy()\ntrack_id = int(tracks_mp.max())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-09T00:16:15.955624Z","iopub.execute_input":"2026-07-09T00:16:15.956667Z","iopub.status.idle":"2026-07-09T00:29:59.742312Z","shell.execute_reply.started":"2026-07-09T00:16:15.956622Z","shell.execute_reply":"2026-07-09T00:29:59.741411Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Con esta funcion se asigna una puntuación de calidad a cada trayectoria reconstruida \n# analizando la coherencia de las conexiones entre sus hits.\ndef get_track_score(tracks_all, n=4):\n    scores = np.zeros(len(tracks_all),dtype=np.float32) # Creamos un array de ceros con la longitud del array de tracks\n    # Recorremos cada trayectoria \n    for i, path in enumerate(tracks_all):\n        count = len(path) # contamos cuantos hits hay\n        # Solo toma en cuenta aquellos con más de un hit\n        if count > 1:\n            tp = 0 # inicializamos los true positives que contarán como trayectoria correcta\n            fp = 0 # inicializamos los false positives que contarán como trayectoria incorrecta\n            \n            # para cada hit\n            for p in path:\n                tp += np.sum(np.isin(tracks_all[p],path,assume_unique=True)) # Cuantos de los hits pertenecen a la trayectoria\n                fp += np.sum(np.isin(tracks_all[p],path,assume_unique=True,invert=True)) # Cuantos no pertenecen a la trayectoria\n            scores[i] = (tp- fp*n- count) / count / (count-1) # Cálculo de la coherencia\n        else:\n            scores[i] = -np.inf\n    return scores\n\n# Devuelve información sobre esas trayectorias, la más destacable el score oficial\ndef evaluate_tracks(tracks, truth):\n    submission = pd.DataFrame({'hit_id':truth.hit_id.values,'track_id':tracks.astype(int)})\n    score = score_event(truth,submission)\n    track_id = tracks.max()\n    print(\"RESULTADOS\")\n    print(f\"Score = {score*100:.2f}%\")\n    print(f\"Tracks reconstruidas: {tracks.max()}\")\n    print(f\"Hits asignados: {np.sum(tracks > 0)}\")\n    print(f\"Hits sin asignar: {np.sum(tracks == 0)}\")\n\n# Función para extender trayectorias. Se intentan alargar.\n# Es muy similar a get_path2 en su primera parte\ndef extend_path(path,mask,thr,last=False):\n    a = 0\n    \n    for p in path[:-1]: # Recorremos todos los hits de la trayectoria salvo el ultimo\n        c = get_predict2(p) # probabilidades\n        \n        if last == False:\n            mask = (c > thr) * mask # Se eliminan los candidatos que no superen el umbral\n        mask[p] = 0 # Se elimina ese hit de la trayectoria\n        cand = np.where( c > thr)[0]\n        if len(cand) > 0:\n            mask[cand[np.isin(module_id[cand],module_id[path])]] = 0 # Se eliminan los del mismo módulo\n        a = (c + a) * mask\n    \n    # Aquí se extiende\n    while True:\n        c = get_predict2(path[-1]) # solo estudiamos las conexiones del último hit\n        \n        if last == False:\n            mask = (c > thr) * mask\n        mask[path[-1]] = 0\n        cand = np.where(c > thr)[0]\n        \n        if len(cand) > 0:\n            mask[cand[np.isin(module_id[cand],module_id[path])]] = 0\n        a = (c + a) * mask \n        \n        # Si el candidato es malo nos salimos y sino se añade al camino\n        if a.max() < thr * len(path):\n            break  \n        path.append(int(a.argmax()))\n        \n        if last:\n            break\n            \n    return path\n    \n\nscores = get_track_score(tracks_all) # Se obtiene la 'calidad' de la trayectoria ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-09T00:31:47.278992Z","iopub.execute_input":"2026-07-09T00:31:47.280090Z","iopub.status.idle":"2026-07-09T00:32:08.902509Z","shell.execute_reply.started":"2026-07-09T00:31:47.280017Z","shell.execute_reply":"2026-07-09T00:32:08.901307Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# MULTIPASO\ntracks_mp = np.zeros(len(hits),dtype=np.int32) # array de ceros\ntrack_id = 0 # inicializador de indice a cero\nidx = np.argsort(scores)[::-1] # Se ordenan de menor a mayor calidad\n# TRAYECTORIAS LARGAS\nfor hit in idx:\n    path = np.array(tracks_all[hit]) # obtenemos el camino para cada hit\n    path = path[np.where(tracks_mp[path] == 0)[0]] # eliminamos hits que ya estan en otras trayectorias\n\n    # Nos quedamos solo con aquellos con más de 6 hits\n    if len(path) > 6:\n        track_id += 1\n        tracks_mp[path] = track_id\nprint(\"TRAYECTORIAS LARGAS\")\nevaluate_tracks(tracks_mp, truth)\n\n\n# Repetimos el proceso extendiendo las trayectorias mayores a 3 hits\nfor hit in idx:\n    path = np.array(tracks_all[hit])\n    path = path[np.where(tracks_mp[path] == 0)[0]]\n    if len(path) > 3:\n        path = extend_path(path.tolist(),1 * (tracks_mp == 0),0.6) # Añadimos aquellos hits que no pertenezcana ninguna \n        track_id += 1\n        tracks_mp[path] = track_id\n\nprint(\"TRAYECTORIAS NUEVAS\")\nevaluate_tracks(tracks_mp, truth)\n\n\n#TRAYECTORIAS CORTAS\n\n# Ahora buscamos extender aquellas con 2 o más hits, solo se conservarán si presentan más de 2 hits\nfor hit in idx:\n    path = np.array(tracks_all[hit])\n    path = path[np.where(tracks_mp[path] == 0)[0]]\n\n    if len(path) > 1:\n        path = extend_path(path.tolist(),1 * (tracks_mp == 0),0.5)\n\n    if len(path) > 2:\n        track_id += 1\n        tracks_mp[path] = track_id\n\nprint(\"TRAYECTORIAS CORTAS\")\nevaluate_tracks(tracks_mp, truth)\n\n# Finalmente se recorren todas las trayectorias de la reconstrucción \nfor tid in range(1,int(tracks_mp.max()) + 1):\n    path = np.where(tracks_mp == tid)[0]  # Se recuperan todos los hits\n\n    # Finalmente la ultima comprobación es intentar alargar aquellas cuyo numero de hits sea par\n    if len(path) % 2 == 0:\n        path = extend_path(path.tolist(),1 * (tracks_mp == 0),0.5,True)\n        tracks_mp[path] = tid\n\nprint(\"FINAL\")\nevaluate_tracks(tracks_mp, truth)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-09T00:32:10.051113Z","iopub.execute_input":"2026-07-09T00:32:10.051911Z","iopub.status.idle":"2026-07-09T00:32:20.652581Z","shell.execute_reply.started":"2026-07-09T00:32:10.051874Z","shell.execute_reply":"2026-07-09T00:32:20.651631Z"}},"outputs":[],"execution_count":null}]}