From 299191024af12ec8b6abd6c505df152f4c1da0a7 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Wed, 3 Jun 2026 17:35:35 +0200 Subject: [PATCH 1/5] Minor. --- src/hexsample/clustering.py | 7 ++++++- src/hexsample/tasks.py | 13 ++++++++++--- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py index 1ee896f..8af34f7 100644 --- a/src/hexsample/clustering.py +++ b/src/hexsample/clustering.py @@ -297,12 +297,17 @@ def run(self, event) -> Optional[Cluster]: # and applying pedestal and gain correction. pha = (event.pha[adc_channel_order] - pedestal(col, row)) / gain(col, row) elif isinstance(event, DigiEventRectangular): + rows, cols = np.mgrid[event.roi.min_row:event.roi.max_row + 1, + event.roi.min_col:event.roi.max_col + 1] + cols = cols.flatten() + rows = rows.flatten() + event.pha = event.pha - pedestal(cols.ravel(), rows.ravel()).reshape(event.roi.shape()) seed_coords = event.highest_pixel() if self.readout.is_at_border(*seed_coords): return None neigh_coords = self.readout.neighbors(*seed_coords) col, row = np.vstack((seed_coords, neigh_coords)).T - pha = (event(col, row) - pedestal(col, row)) / gain(col, row) + pha = event(col, row) / gain(col, row) else: raise RuntimeError(f"Unsupported event type {type(event)} for clustering") # Zero suppressing the event (whatever the readout type)... diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 5f631c8..b4f5ecb 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -997,16 +997,23 @@ def calibview( f"lower_quantile={lower_quantile}, upper_quantile={upper_quantile}" ) logger.info(f"Number of calibrated pixels after quality cuts: {np.sum(mask)}") + mask = mask & (matrix.values >= lower_bound) & (matrix.values <= upper_bound) # Plot the values matrix. plt.figure(f"Calibrated matrix: {matrix.metadata['file_name']}") plt.imshow(matrix.values, origin="upper", vmin=lower_bound, vmax=upper_bound) plt.xlabel("Column") plt.ylabel("Row") plt.colorbar(label=unit) + # Plot the entries matrix. + plt.figure(f"Entries matrix: {matrix.metadata['file_name']}") + plt.imshow(matrix.entries, origin="upper") + plt.xlabel("Column") + plt.ylabel("Row") + plt.colorbar(label="Entries") # Plot the distribution of the calibrated values. vals = matrix.values.flatten()[mask.flatten()] edges = np.linspace(lower_bound, upper_bound, 100) - vals_hist = Histogram1d(edges, label="Distribution", xlabel=unit).fill(vals) + vals_hist = Histogram1d(edges, label="Values distribution", xlabel=unit).fill(vals) plt.figure("Distribution of calibrated values") vals_hist.plot(statistics=True) plt.legend() @@ -1051,10 +1058,10 @@ def calibview( plt.xlabel(f"Calibrated values [{unit}]") plt.ylabel(f"Monte Carlo truth values [{mc_unit}]") # Plot the residuals distribution. - residuals = (vals - mc_vals[mask.flatten()]) / mc_vals[mask.flatten()] + residuals = vals - mc_vals[mask.flatten()] residual_edges = np.linspace(np.nanmin(residuals), np.nanmax(residuals), 100) residual_hist = Histogram1d( - residual_edges, label="Residuals", xlabel="Relative Residual" + residual_edges, label="Residuals distribution", xlabel=f"Residuals [{unit}]" ).fill(residuals) plt.figure("Relative residuals distribution") residual_hist.plot(statistics=True) From 029c55e3bae3e7d2730e3c40d99e0de4816e6a10 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Wed, 3 Jun 2026 18:16:37 +0200 Subject: [PATCH 2/5] Test. --- src/hexsample/likelihood.py | 195 +++++++++++++++++++++++++++++------- src/hexsample/position.py | 150 +++++++++++++++++---------- src/hexsample/tasks.py | 2 +- 3 files changed, 260 insertions(+), 87 deletions(-) diff --git a/src/hexsample/likelihood.py b/src/hexsample/likelihood.py index fd70d29..53c614f 100644 --- a/src/hexsample/likelihood.py +++ b/src/hexsample/likelihood.py @@ -186,51 +186,178 @@ def weighted_pha(pha: np.ndarray, f_interp: np.ndarray, inv_sigma2: np.ndarray) return max(.0, sum_qf / (sum_f2 + 1e-12)) -@njit -def nll_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, - bin_size: float, noise: np.ndarray) -> float: - """Compute the negative log-likelihood for a given position (x, y). +# @njit +# def nll_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, +# bin_size: float, noise: np.ndarray) -> float: +# """Compute the negative log-likelihood for a given position (x, y). + +# The model is based on the Gaussian diffusion of the charge cloud, and uses the precomputed +# charge fractions in each pixel from the f map. The summed pha is profiled out to reduce the +# dimensionality of the optimization. +# """ +# # Calculate the bin indices and fractional coordinates for the interpolation +# ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) +# # Interpolate the charge fractions for the 7 pixels in the cluster +# f_interp = interpolation(f, ix0, iy0, wx, wy) +# # Calculate the inverse of the noise variance for each pixel +# inv_sigma2 = 1.0 / (noise**2) +# # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) +# total_pha = weighted_pha(pha, f_interp, inv_sigma2) +# # Now compute the NLL using the optimal energy +# nll = 0.0 +# for i in range(7): +# mu = f_interp[i] * total_pha +# res = pha[i] - mu +# nll += 0.5 * (res**2 * inv_sigma2[i] + LOG2PI) +# return nll + - The model is based on the Gaussian diffusion of the charge cloud, and uses the precomputed - charge fractions in each pixel from the f map. The summed pha is profiled out to reduce the - dimensionality of the optimization. +# @njit +# def nll_grad_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, +# bin_size: float, noise: np.ndarray) -> np.ndarray: +# """Compute the gradient of the negative log-likelihood with respect to the free parameters. +# """ +# # Calculate the bin indices and fractional coordinates for the interpolation +# ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) +# # Interpolate the charge fractions and their derivatives for the 7 pixels in the cluster +# f_interp, df_dx, df_dy = interpolation_derivatives(f, ix0, iy0, wx, wy, bin_size) +# # Calculate the inverse of the noise variance for each pixel +# inv_sigma2 = 1.0 / (noise**2) +# # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) +# total_pha = weighted_pha(pha, f_interp, inv_sigma2) +# # Now compute the gradient using the optimal energy +# gnll_x = 0.0 +# gnll_y = 0.0 +# for i in range(7): +# mu = f_interp[i] * total_pha +# d_loss_dmu = -(pha[i] - mu) * inv_sigma2[i] +# gnll_x += d_loss_dmu * total_pha * df_dx[i] +# gnll_y += d_loss_dmu * total_pha * df_dy[i] +# return np.array([gnll_x, gnll_y]) + + +@njit +def nll_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, + xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> float: + """Compute the exact multivariate Gauss-Multinomial negative log-likelihood. + Includes both electronic noise and physical charge sharing covariance. """ - # Calculate the bin indices and fractional coordinates for the interpolation + # Evitiamo valori fisicamente impossibili per la carica totale nell'ottimizzatore + if total_pha <= 0.0: + return 1e12 + + # 1. Coordinate e interpolazione geometrica ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) - # Interpolate the charge fractions for the 7 pixels in the cluster f_interp = interpolation(f, ix0, iy0, wx, wy) - # Calculate the inverse of the noise variance for each pixel - inv_sigma2 = 1.0 / (noise**2) - # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) - total_pha = weighted_pha(pha, f_interp, inv_sigma2) - # Now compute the NLL using the optimal energy - nll = 0.0 + + # 2. Costruzione analitica dei vettori di supporto per Sherman-Morrison + # D_ii = sigma_e^2 + Q * f_i + D = np.zeros(7) + inv_D = np.zeros(7) + for i in range(7): + D[i] = (noise[i] ** 2) + total_pha * f_interp[i] + inv_D[i] = 1.0 / D[i] + + # Calcolo del fattore di accoppiamento multinomiale: sum( f_k^2 / D_kk ) + somma_fk2_D = 0.0 + for k in range(7): + somma_fk2_D += (f_interp[k] ** 2) * inv_D[k] + + # Il fattore correttivo dell'inversa (denominatore di Sherman-Morrison) + # Nota: aggiungiamo un segno meno coerente con la covarianza negativa dello sharing + g = 1.0 / (1.0 - total_pha * somma_fk2_D + 1e-12) + + # 3. Calcolo del Log-Determinante di V per il termine di normalizzazione della likelihood + prod_D = 0.0 + for k in range(7): + prod_D += math.log(D[k]) + # det(V) = prod(D) * (1 - Q * sum(f^2/D)) + log_det_V = prod_D + math.log(abs(1.0 - total_pha * somma_fk2_D) + 1e-12) + + # 4. Calcolo del vettore dei residui (Misurato - Atteso) + res = np.zeros(7) + for i in range(7): + res[i] = pha[i] - (total_pha * f_interp[i]) + + # 5. Moltiplicazione del vettore residui per la matrice inversa analitica: res^T * V^-1 * res + # Sviluppando l'algebra di Sherman-Morrison, non serve allocare la matrice 7x7! + termine_diagonale_res = 0.0 + termine_incrociato_res = 0.0 + for i in range(7): - mu = f_interp[i] * total_pha - res = pha[i] - mu - nll += 0.5 * (res**2 * inv_sigma2[i] + LOG2PI) + termine_diagonale_res += (res[i] ** 2) * inv_D[i] + termine_incrociato_res += res[i] * f_interp[i] * inv_D[i] + + chi2_multivariato = termine_diagonale_res + g * total_pha * (termine_incrociato_res ** 2) + + # 6. Somma finale della NLL + nll = 0.5 * chi2_multivariato + 0.5 * log_det_V + (7.2 / 2.0) * LOG2PI + return nll - @njit -def nll_grad_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, - bin_size: float, noise: np.ndarray) -> np.ndarray: - """Compute the gradient of the negative log-likelihood with respect to the free parameters. +def nll_grad_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, + xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> np.ndarray: + """Compute the complete analytical gradient of the exact multivariate NLL + with respect to x, y, and Q simultaneously (Returns array of size 3). """ - # Calculate the bin indices and fractional coordinates for the interpolation + if total_pha <= 0.0: + return np.array([0.0, 0.0, 0.0]) + + # 1. Una sola chiamata geometrica per frazioni e derivate spaziali ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) - # Interpolate the charge fractions and their derivatives for the 7 pixels in the cluster f_interp, df_dx, df_dy = interpolation_derivatives(f, ix0, iy0, wx, wy, bin_size) - # Calculate the inverse of the noise variance for each pixel - inv_sigma2 = 1.0 / (noise**2) - # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) - total_pha = weighted_pha(pha, f_interp, inv_sigma2) - # Now compute the gradient using the optimal energy + + # 2. Ricostruzione del vettore D e del nucleo di Sherman-Morrison (Condiviso) + D = np.zeros(7) + inv_D = np.zeros(7) + for i in range(7): + D[i] = (noise[i] ** 2) + total_pha * f_interp[i] + inv_D[i] = 1.0 / D[i] + + somma_fk2_D = 0.0 + for k in range(7): + somma_fk2_D += (f_interp[k] ** 2) * inv_D[k] + + g = 1.0 / (1.0 - total_pha * somma_fk2_D + 1e-12) + + # 3. Calcolo dei residui e del termine incrociato (Condiviso) + res = np.zeros(7) + termine_incrociato_res = 0.0 + for i in range(7): + res[i] = pha[i] - (total_pha * f_interp[i]) + termine_incrociato_res += res[i] * f_interp[i] * inv_D[i] + + # ========================================================================= + # PARTE Spaziale (x, y) + # ========================================================================= + alpha = g * total_pha * termine_incrociato_res + vettore_w = np.zeros(7) + for i in range(7): + vettore_w[i] = inv_D[i] * (res[i] + alpha * f_interp[i]) + + dNLL_df = np.zeros(7) + for i in range(7): + dNLL_df[i] = total_pha * ( + 0.5 * inv_D[i] + - vettore_w[i] + - 0.5 * (vettore_w[i] ** 2) + + g * vettore_w[i] * f_interp[i] * (vettore_w[i] - inv_D[i] * res[i]) + ) + gnll_x = 0.0 gnll_y = 0.0 for i in range(7): - mu = f_interp[i] * total_pha - d_loss_dmu = -(pha[i] - mu) * inv_sigma2[i] - gnll_x += d_loss_dmu * total_pha * df_dx[i] - gnll_y += d_loss_dmu * total_pha * df_dy[i] - return np.array([gnll_x, gnll_y]) + gnll_x += dNLL_df[i] * df_dx[i] + gnll_y += dNLL_df[i] * df_dy[i] + + # ========================================================================= + # PARTE Energetica (Q) - Calcolata sfruttando gli stessi loop + # ========================================================================= + gnll_q = 0.0 + for i in range(7): + # Derivata rispetto a Q coerente col peso statistico del rumore dei pixel + gnll_q -= (res[i] * f_interp[i]) / (noise[i] ** 2) + + # Ritorna l'array nativo a 3 dimensioni [grad_x, grad_y, grad_q] + return np.array([gnll_x, gnll_y, gnll_q]) \ No newline at end of file diff --git a/src/hexsample/position.py b/src/hexsample/position.py index 370dbb0..bf45522 100644 --- a/src/hexsample/position.py +++ b/src/hexsample/position.py @@ -269,6 +269,95 @@ def eta_3pix( return dx, dy +# def mle( +# pha: np.ndarray, +# noise: np.ndarray, +# f: np.ndarray, +# bin_size: float, +# xlims: Tuple[float, float], +# ylims: Tuple[float, float], +# p0: Tuple[float, float] = (0.0, 0.0), +# ) -> Minuit: +# """Perform maximum likelihood estimation of the incident position of the +# photon, given the observed pha in the 7 pixels of the cluster. + +# The likelihood used for the fit is based on the Gaussian diffusion of the +# charge cloud around the incident position, and uses the precomputed charge +# fractions in each pixel to evaluate the likelihood for a given position. + +# To speed up the computation, the negative log-likelihood and its gradient +# are implemented in the likelihood.py module and decorated with numba.njit. + +# Arguments +# --------- +# pha : np.ndarray +# The measured pha in the 7 pixels of the cluster, ordered according to +# the convention defined in calibration.py. + +# noise : np.ndarray +# The array of shape (7,) containing the equalized noise standard deviation +# for each pixel. + +# f : np.ndarray +# The array of shape (7, nx, ny) containing the precomputed charge fractions +# in each pixel as a function of the incident position. + +# bin_size : float +# The size of the bins in the f array, expressed in units of the pixel +# pitch. + +# xlims : Tuple[float, float] +# The limits for the x coordinate of the f array, expressed in +# units of the pixel pitch. + +# ylims : Tuple[float, float] +# The limits for the y coordinate of the f array, expressed in +# units of the pixel pitch. + +# p0 : Tuple[float, float], optional +# The initial guess for the (x, y) position of the photon, expressed in +# units of the pixel pitch. A reasonable initial guess can be the centroid +# of the cluster. Default is the center of the pixel (0.0, 0.0). + +# Returns +# ------- +# m : Minuit +# The minimizer object containing all the information about the fit. +# """ +# # Unpack the grid limits on the x and y axes. +# xmin, xmax = xlims +# ymin, ymax = ylims +# # Define the objective functions for the optimization, which are the +# # negative log-likelihood... +# def nll(x: float, y: float) -> float: +# return nll_numba(x, y, pha, f, xmin, ymin, bin_size, noise) +# # ... and its gradient. +# def nll_grad(x: float, y: float) -> Tuple[float, float]: +# return nll_grad_numba(x, y, pha, f, xmin, ymin, bin_size, noise) +# # Assign a name to the free parameters. +# parnames = ["x", "y"] +# # Initialize the minimizer. +# m = Minuit(nll, *p0, grad=nll_grad, name=parnames) +# # Set the limits for the free parameters. +# m.limits["x"] = (xmin, xmax) +# m.limits["y"] = (ymin, ymax) +# # Set the initial step sizes for the minimizer. We choose half the bin +# # size as default value. This value is automatically adjusted by the +# # minimizer during the fit, so it is not critical to choose a very +# # precise value. +# m.errors["x"] = bin_size / 2 +# m.errors["y"] = bin_size / 2 +# # Define the strategy for the minimizer. We use the higher strategy to +# # avoid numerical problems with the hessian. For an explanation of the +# # strategy levels, please refer to the iminuit documentation. +# m.strategy = 2 +# # Run the minimization. +# m.migrad() +# # Return the minimizer object, which contains the best-fit values and +# # their uncertainties. It's better to return the whole object to allow +# # the caller to inspect the fit results and diagnostics. +# return m + def mle( pha: np.ndarray, noise: np.ndarray, @@ -278,69 +367,26 @@ def mle( ylims: Tuple[float, float], p0: Tuple[float, float] = (0.0, 0.0), ) -> Minuit: - """Perform maximum likelihood estimation of the incident position of the - photon, given the observed pha in the 7 pixels of the cluster. - - The likelihood used for the fit is based on the Gaussian diffusion of the - charge cloud around the incident position, and uses the precomputed charge - fractions in each pixel to evaluate the likelihood for a given position. - - To speed up the computation, the negative log-likelihood and its gradient - are implemented in the likelihood.py module and decorated with numba.njit. - - Arguments - --------- - pha : np.ndarray - The measured pha in the 7 pixels of the cluster, ordered according to - the convention defined in calibration.py. - - noise : np.ndarray - The array of shape (7,) containing the equalized noise standard deviation - for each pixel. - - f : np.ndarray - The array of shape (7, nx, ny) containing the precomputed charge fractions - in each pixel as a function of the incident position. - - bin_size : float - The size of the bins in the f array, expressed in units of the pixel - pitch. - - xlims : Tuple[float, float] - The limits for the x coordinate of the f array, expressed in - units of the pixel pitch. - - ylims : Tuple[float, float] - The limits for the y coordinate of the f array, expressed in - units of the pixel pitch. - - p0 : Tuple[float, float], optional - The initial guess for the (x, y) position of the photon, expressed in - units of the pixel pitch. A reasonable initial guess can be the centroid - of the cluster. Default is the center of the pixel (0.0, 0.0). - - Returns - ------- - m : Minuit - The minimizer object containing all the information about the fit. - """ + p0 = (*p0, np.sum(pha)) # Unpack the grid limits on the x and y axes. xmin, xmax = xlims ymin, ymax = ylims # Define the objective functions for the optimization, which are the # negative log-likelihood... - def nll(x: float, y: float) -> float: - return nll_numba(x, y, pha, f, xmin, ymin, bin_size, noise) + def nll(x: float, y: float, q: float) -> float: + return nll_numba(x, y, q, pha, f, xmin, ymin, bin_size, noise) # ... and its gradient. - def nll_grad(x: float, y: float) -> Tuple[float, float]: - return nll_grad_numba(x, y, pha, f, xmin, ymin, bin_size, noise) + def nll_grad(x: float, y: float, q: float) -> Tuple[float, float]: + return nll_grad_numba(x, y, q, pha, f, xmin, ymin, bin_size, noise) # Assign a name to the free parameters. - parnames = ["x", "y"] + parnames = ["x", "y", "q"] # Initialize the minimizer. m = Minuit(nll, *p0, grad=nll_grad, name=parnames) # Set the limits for the free parameters. m.limits["x"] = (xmin, xmax) m.limits["y"] = (ymin, ymax) + m.limits["q"] = (np.sum(pha) * 0.5, np.sum(pha) * 1.5) + # m.fixed["q"] = True # Set the initial step sizes for the minimizer. We choose half the bin # size as default value. This value is automatically adjusted by the # minimizer during the fit, so it is not critical to choose a very @@ -356,4 +402,4 @@ def nll_grad(x: float, y: float) -> Tuple[float, float]: # Return the minimizer object, which contains the best-fit values and # their uncertainties. It's better to return the whole object to allow # the caller to inspect the fit results and diagnostics. - return m + return m \ No newline at end of file diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index b4f5ecb..a3c8f05 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -300,7 +300,7 @@ def reconstruct( equalization_matrix=equalization_matrix ) if pos_recon_algorithm == "mle": - clustering = ClusteringHex(readout, 0, pos_recon_algorithm, recon_pars) + clustering = ClusteringHex(readout, -5, pos_recon_algorithm, recon_pars) num_neighbors = 6 else: clustering = ClusteringNN( From ed2ebd6a6ffb288e4f2b3f00d12ab7ee39e58978 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Fri, 19 Jun 2026 09:45:21 +0200 Subject: [PATCH 3/5] Minor. --- src/hexsample/position.py | 1 - src/hexsample/tasks.py | 2 +- tests/test_likelihood.py | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/hexsample/position.py b/src/hexsample/position.py index bf45522..ceb7994 100644 --- a/src/hexsample/position.py +++ b/src/hexsample/position.py @@ -386,7 +386,6 @@ def nll_grad(x: float, y: float, q: float) -> Tuple[float, float]: m.limits["x"] = (xmin, xmax) m.limits["y"] = (ymin, ymax) m.limits["q"] = (np.sum(pha) * 0.5, np.sum(pha) * 1.5) - # m.fixed["q"] = True # Set the initial step sizes for the minimizer. We choose half the bin # size as default value. This value is automatically adjusted by the # minimizer during the fit, so it is not critical to choose a very diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index a3c8f05..68cf1b0 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -416,7 +416,7 @@ def calibrate_position( readout = create_readout(readout_mode, header, *readout_args) # To correctly analyze every type of event, we need a zero suppression threshold # of 0, because the calibration should be performed on zero-noise simulations. - clustering = ClusteringHex(readout, zero_sup_threshold=0.) + clustering = ClusteringHex(readout, zero_sup_threshold=-5.) clustering_nn = ClusteringNN(readout, zero_sup_threshold=zero_sup_threshold, num_neighbors=6, pos_recon_algorithm="centroid") # Initialize the position calibrator and run the event loop. diff --git a/tests/test_likelihood.py b/tests/test_likelihood.py index c193b3f..6618f0b 100644 --- a/tests/test_likelihood.py +++ b/tests/test_likelihood.py @@ -44,7 +44,7 @@ def test_nll_numba(): nll = np.zeros((len(x), len(y))) for i_x, _x in enumerate(x): for i_y, _y in enumerate(y): - nll[i_x, i_y] = nll_numba(_x, _y, pha, f, xbin0, ybin0, bin_size, sigma) + nll[i_x, i_y] = nll_numba(_x, _y, np.sum(pha),pha, f, xbin0, ybin0, bin_size, sigma) plt.figure("test_negative_log_likelihood") plt.imshow(nll.T, extent=(x[0], x[-1], y[0], y[-1]), origin="lower") plt.colorbar(label="Negative log-likelihood") From 0d66585e8f69f7a27c95a8a23988b19324ab9b13 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Sun, 26 Jul 2026 12:11:42 +0200 Subject: [PATCH 4/5] MLE fixed. --- src/hexsample/likelihood.py | 227 ++++++++++++------------------------ src/hexsample/position.py | 142 ++++++++-------------- 2 files changed, 127 insertions(+), 242 deletions(-) diff --git a/src/hexsample/likelihood.py b/src/hexsample/likelihood.py index 53c614f..524cab0 100644 --- a/src/hexsample/likelihood.py +++ b/src/hexsample/likelihood.py @@ -186,178 +186,105 @@ def weighted_pha(pha: np.ndarray, f_interp: np.ndarray, inv_sigma2: np.ndarray) return max(.0, sum_qf / (sum_f2 + 1e-12)) -# @njit -# def nll_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, -# bin_size: float, noise: np.ndarray) -> float: -# """Compute the negative log-likelihood for a given position (x, y). - -# The model is based on the Gaussian diffusion of the charge cloud, and uses the precomputed -# charge fractions in each pixel from the f map. The summed pha is profiled out to reduce the -# dimensionality of the optimization. -# """ -# # Calculate the bin indices and fractional coordinates for the interpolation -# ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) -# # Interpolate the charge fractions for the 7 pixels in the cluster -# f_interp = interpolation(f, ix0, iy0, wx, wy) -# # Calculate the inverse of the noise variance for each pixel -# inv_sigma2 = 1.0 / (noise**2) -# # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) -# total_pha = weighted_pha(pha, f_interp, inv_sigma2) -# # Now compute the NLL using the optimal energy -# nll = 0.0 -# for i in range(7): -# mu = f_interp[i] * total_pha -# res = pha[i] - mu -# nll += 0.5 * (res**2 * inv_sigma2[i] + LOG2PI) -# return nll - +@njit +def sherman_morrison_inverse(res: np.ndarray, f_interp: np.ndarray, total_pha: float, + noise: np.ndarray): + """Compute the chi2 (res^T * V^-1 * res) and log-determinant of the covariance + matrix V using the Sherman-Morrison formula and the matrix determinant lemma. -# @njit -# def nll_grad_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, -# bin_size: float, noise: np.ndarray) -> np.ndarray: -# """Compute the gradient of the negative log-likelihood with respect to the free parameters. -# """ -# # Calculate the bin indices and fractional coordinates for the interpolation -# ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) -# # Interpolate the charge fractions and their derivatives for the 7 pixels in the cluster -# f_interp, df_dx, df_dy = interpolation_derivatives(f, ix0, iy0, wx, wy, bin_size) -# # Calculate the inverse of the noise variance for each pixel -# inv_sigma2 = 1.0 / (noise**2) -# # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) -# total_pha = weighted_pha(pha, f_interp, inv_sigma2) -# # Now compute the gradient using the optimal energy -# gnll_x = 0.0 -# gnll_y = 0.0 -# for i in range(7): -# mu = f_interp[i] * total_pha -# d_loss_dmu = -(pha[i] - mu) * inv_sigma2[i] -# gnll_x += d_loss_dmu * total_pha * df_dx[i] -# gnll_y += d_loss_dmu * total_pha * df_dy[i] -# return np.array([gnll_x, gnll_y]) + This avoids the explicit inversion of the 7x7 matrix and allows to save space and time. + """ + # First we compute the diagonal elements of D, given by the square sum of electronic + # noise and statistical fluctuations, and their inverses. + d_matrix = np.zeros(7) + inv_d_matrix = np.zeros(7) + for i in range(7): + d_matrix[i] = (noise[i] ** 2) + total_pha * f_interp[i] + inv_d_matrix[i] = 1.0 / d_matrix[i] + # Then we calculate the term at the denominator. + sum_fk2_d = 0.0 + for k in range(7): + sum_fk2_d += (f_interp[k] ** 2) * inv_d_matrix[k] + # The final denominator term. + g = 1.0 / (1.0 - total_pha * sum_fk2_d + 1e-12) + # Calculate the log determinant of V using the matrix determinant lemma. + prod_d = 0.0 + for k in range(7): + prod_d += math.log(d_matrix[k]) + log_det_v = prod_d + math.log(abs(1.0 - total_pha * sum_fk2_d) + 1e-12) + # Calculate the quadratic form res^T * V^-1 * res using Sherman-Morrison algebra. + diagional_term = 0.0 + cross_term = 0.0 + for i in range(7): + diagional_term += (res[i] ** 2) * inv_d_matrix[i] + cross_term += res[i] * f_interp[i] * inv_d_matrix[i] + # Final chi2 value. + chi2 = diagional_term + g * total_pha * (cross_term ** 2) + # Calculate the vector w = V^-1 * res for the gradient calculation. + alpha = g * total_pha * cross_term + vettore_w = np.zeros(7) + for i in range(7): + vettore_w[i] = inv_d_matrix[i] * (res[i] + alpha * f_interp[i]) + dnll_df = np.zeros(7) + for i in range(7): + dnll_df[i] = total_pha * ( + 0.5 * inv_d_matrix[i] + - vettore_w[i] + - 0.5 * (vettore_w[i] ** 2) + + g * vettore_w[i] * f_interp[i] * (vettore_w[i] - inv_d_matrix[i] * res[i]) + ) + return chi2, log_det_v, dnll_df @njit -def nll_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, - xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> float: - """Compute the exact multivariate Gauss-Multinomial negative log-likelihood. - Includes both electronic noise and physical charge sharing covariance. +def nll_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, + xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> float: + """Compute the negative log-likelihood for a given position (x, y). + + The model is based on the Gaussian diffusion of the charge cloud, and uses the precomputed + charge fractions in each pixel from the f map. """ - # Evitiamo valori fisicamente impossibili per la carica totale nell'ottimizzatore - if total_pha <= 0.0: - return 1e12 - - # 1. Coordinate e interpolazione geometrica + # Calculate the bin indices and fractional coordinates for the interpolation ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) + # Interpolate the charge fractions for the 7 pixels in the cluster f_interp = interpolation(f, ix0, iy0, wx, wy) - - # 2. Costruzione analitica dei vettori di supporto per Sherman-Morrison - # D_ii = sigma_e^2 + Q * f_i - D = np.zeros(7) - inv_D = np.zeros(7) - for i in range(7): - D[i] = (noise[i] ** 2) + total_pha * f_interp[i] - inv_D[i] = 1.0 / D[i] - - # Calcolo del fattore di accoppiamento multinomiale: sum( f_k^2 / D_kk ) - somma_fk2_D = 0.0 - for k in range(7): - somma_fk2_D += (f_interp[k] ** 2) * inv_D[k] - - # Il fattore correttivo dell'inversa (denominatore di Sherman-Morrison) - # Nota: aggiungiamo un segno meno coerente con la covarianza negativa dello sharing - g = 1.0 / (1.0 - total_pha * somma_fk2_D + 1e-12) - - # 3. Calcolo del Log-Determinante di V per il termine di normalizzazione della likelihood - prod_D = 0.0 - for k in range(7): - prod_D += math.log(D[k]) - # det(V) = prod(D) * (1 - Q * sum(f^2/D)) - log_det_V = prod_D + math.log(abs(1.0 - total_pha * somma_fk2_D) + 1e-12) - - # 4. Calcolo del vettore dei residui (Misurato - Atteso) + # We need to invert the covariance matrix V. We use the Sherman-Morrison formula + # to avoid inverting a full 7x7 matrix. This is possible because V is a diagonal + # matrix D plus a rank 1 matrix. res = np.zeros(7) for i in range(7): res[i] = pha[i] - (total_pha * f_interp[i]) - - # 5. Moltiplicazione del vettore residui per la matrice inversa analitica: res^T * V^-1 * res - # Sviluppando l'algebra di Sherman-Morrison, non serve allocare la matrice 7x7! - termine_diagonale_res = 0.0 - termine_incrociato_res = 0.0 - - for i in range(7): - termine_diagonale_res += (res[i] ** 2) * inv_D[i] - termine_incrociato_res += res[i] * f_interp[i] * inv_D[i] - - chi2_multivariato = termine_diagonale_res + g * total_pha * (termine_incrociato_res ** 2) - - # 6. Somma finale della NLL - nll = 0.5 * chi2_multivariato + 0.5 * log_det_V + (7.2 / 2.0) * LOG2PI - + chi2, log_det_v, _ = sherman_morrison_inverse(res, f_interp, total_pha, noise) + # Finally we compute the negative log-likelihood using the chi2 and the log-determinant of V. + nll = 0.5 * chi2 + 0.5 * log_det_v return nll + @njit -def nll_grad_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, - xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> np.ndarray: - """Compute the complete analytical gradient of the exact multivariate NLL - with respect to x, y, and Q simultaneously (Returns array of size 3). +def nll_grad_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, + xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> np.ndarray: + """Compute the complete analytical gradient of the exact multivariate NLL with + respect to x, y, and Q simultaneously (Returns array of size 3). """ - if total_pha <= 0.0: - return np.array([0.0, 0.0, 0.0]) - - # 1. Una sola chiamata geometrica per frazioni e derivate spaziali + # Calculate the bin indices and fractional coordinates for the interpolation ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) + # Interpolate the charge fractions and their derivatives for the 7 pixels in the cluster f_interp, df_dx, df_dy = interpolation_derivatives(f, ix0, iy0, wx, wy, bin_size) - - # 2. Ricostruzione del vettore D e del nucleo di Sherman-Morrison (Condiviso) - D = np.zeros(7) - inv_D = np.zeros(7) - for i in range(7): - D[i] = (noise[i] ** 2) + total_pha * f_interp[i] - inv_D[i] = 1.0 / D[i] - - somma_fk2_D = 0.0 - for k in range(7): - somma_fk2_D += (f_interp[k] ** 2) * inv_D[k] - - g = 1.0 / (1.0 - total_pha * somma_fk2_D + 1e-12) - - # 3. Calcolo dei residui e del termine incrociato (Condiviso) + # Calculate the residuals for the 7 pixels in the cluster. res = np.zeros(7) - termine_incrociato_res = 0.0 for i in range(7): res[i] = pha[i] - (total_pha * f_interp[i]) - termine_incrociato_res += res[i] * f_interp[i] * inv_D[i] - - # ========================================================================= - # PARTE Spaziale (x, y) - # ========================================================================= - alpha = g * total_pha * termine_incrociato_res - vettore_w = np.zeros(7) - for i in range(7): - vettore_w[i] = inv_D[i] * (res[i] + alpha * f_interp[i]) - - dNLL_df = np.zeros(7) - for i in range(7): - dNLL_df[i] = total_pha * ( - 0.5 * inv_D[i] - - vettore_w[i] - - 0.5 * (vettore_w[i] ** 2) - + g * vettore_w[i] * f_interp[i] * (vettore_w[i] - inv_D[i] * res[i]) - ) - + # We need to invert the covariance matrix V. We use the Sherman-Morrison formula + # to avoid inverting a full 7x7 matrix. + _, _, dnll_df = sherman_morrison_inverse(res, f_interp, total_pha, noise) + # Finally we compute the gradient of the negative log-likelihood with respect to x, y, and Q. gnll_x = 0.0 gnll_y = 0.0 for i in range(7): - gnll_x += dNLL_df[i] * df_dx[i] - gnll_y += dNLL_df[i] * df_dy[i] - - # ========================================================================= - # PARTE Energetica (Q) - Calcolata sfruttando gli stessi loop - # ========================================================================= + gnll_x += dnll_df[i] * df_dx[i] + gnll_y += dnll_df[i] * df_dy[i] gnll_q = 0.0 for i in range(7): - # Derivata rispetto a Q coerente col peso statistico del rumore dei pixel gnll_q -= (res[i] * f_interp[i]) / (noise[i] ** 2) - - # Ritorna l'array nativo a 3 dimensioni [grad_x, grad_y, grad_q] - return np.array([gnll_x, gnll_y, gnll_q]) \ No newline at end of file + # Return the gradient as a numpy array of size 3. + return np.array([gnll_x, gnll_y, gnll_q]) diff --git a/src/hexsample/position.py b/src/hexsample/position.py index ceb7994..4a7a13a 100644 --- a/src/hexsample/position.py +++ b/src/hexsample/position.py @@ -269,95 +269,6 @@ def eta_3pix( return dx, dy -# def mle( -# pha: np.ndarray, -# noise: np.ndarray, -# f: np.ndarray, -# bin_size: float, -# xlims: Tuple[float, float], -# ylims: Tuple[float, float], -# p0: Tuple[float, float] = (0.0, 0.0), -# ) -> Minuit: -# """Perform maximum likelihood estimation of the incident position of the -# photon, given the observed pha in the 7 pixels of the cluster. - -# The likelihood used for the fit is based on the Gaussian diffusion of the -# charge cloud around the incident position, and uses the precomputed charge -# fractions in each pixel to evaluate the likelihood for a given position. - -# To speed up the computation, the negative log-likelihood and its gradient -# are implemented in the likelihood.py module and decorated with numba.njit. - -# Arguments -# --------- -# pha : np.ndarray -# The measured pha in the 7 pixels of the cluster, ordered according to -# the convention defined in calibration.py. - -# noise : np.ndarray -# The array of shape (7,) containing the equalized noise standard deviation -# for each pixel. - -# f : np.ndarray -# The array of shape (7, nx, ny) containing the precomputed charge fractions -# in each pixel as a function of the incident position. - -# bin_size : float -# The size of the bins in the f array, expressed in units of the pixel -# pitch. - -# xlims : Tuple[float, float] -# The limits for the x coordinate of the f array, expressed in -# units of the pixel pitch. - -# ylims : Tuple[float, float] -# The limits for the y coordinate of the f array, expressed in -# units of the pixel pitch. - -# p0 : Tuple[float, float], optional -# The initial guess for the (x, y) position of the photon, expressed in -# units of the pixel pitch. A reasonable initial guess can be the centroid -# of the cluster. Default is the center of the pixel (0.0, 0.0). - -# Returns -# ------- -# m : Minuit -# The minimizer object containing all the information about the fit. -# """ -# # Unpack the grid limits on the x and y axes. -# xmin, xmax = xlims -# ymin, ymax = ylims -# # Define the objective functions for the optimization, which are the -# # negative log-likelihood... -# def nll(x: float, y: float) -> float: -# return nll_numba(x, y, pha, f, xmin, ymin, bin_size, noise) -# # ... and its gradient. -# def nll_grad(x: float, y: float) -> Tuple[float, float]: -# return nll_grad_numba(x, y, pha, f, xmin, ymin, bin_size, noise) -# # Assign a name to the free parameters. -# parnames = ["x", "y"] -# # Initialize the minimizer. -# m = Minuit(nll, *p0, grad=nll_grad, name=parnames) -# # Set the limits for the free parameters. -# m.limits["x"] = (xmin, xmax) -# m.limits["y"] = (ymin, ymax) -# # Set the initial step sizes for the minimizer. We choose half the bin -# # size as default value. This value is automatically adjusted by the -# # minimizer during the fit, so it is not critical to choose a very -# # precise value. -# m.errors["x"] = bin_size / 2 -# m.errors["y"] = bin_size / 2 -# # Define the strategy for the minimizer. We use the higher strategy to -# # avoid numerical problems with the hessian. For an explanation of the -# # strategy levels, please refer to the iminuit documentation. -# m.strategy = 2 -# # Run the minimization. -# m.migrad() -# # Return the minimizer object, which contains the best-fit values and -# # their uncertainties. It's better to return the whole object to allow -# # the caller to inspect the fit results and diagnostics. -# return m - def mle( pha: np.ndarray, noise: np.ndarray, @@ -367,7 +278,52 @@ def mle( ylims: Tuple[float, float], p0: Tuple[float, float] = (0.0, 0.0), ) -> Minuit: - p0 = (*p0, np.sum(pha)) + """Perform maximum likelihood estimation of the incident position of the + photon, given the observed pha in the 7 pixels of the cluster. + + The likelihood used for the fit is based on the Gaussian diffusion of the + charge cloud around the incident position, and uses the precomputed charge + fractions in each pixel to evaluate the likelihood for a given position. + + To speed up the computation, the negative log-likelihood and its gradient + are implemented in the likelihood.py module and decorated with numba.njit. + + Arguments + --------- + pha : np.ndarray + The measured pha in the 7 pixels of the cluster, ordered according to + the convention defined in calibration.py. + + noise : np.ndarray + The array of shape (7,) containing the equalized noise standard deviation + for each pixel. + + f : np.ndarray + The array of shape (7, nx, ny) containing the precomputed charge fractions + in each pixel as a function of the incident position. + + bin_size : float + The size of the bins in the f array, expressed in units of the pixel + pitch. + + xlims : Tuple[float, float] + The limits for the x coordinate of the f array, expressed in + units of the pixel pitch. + + ylims : Tuple[float, float] + The limits for the y coordinate of the f array, expressed in + units of the pixel pitch. + + p0 : Tuple[float, float], optional + The initial guess for the (x, y) position of the photon, expressed in + units of the pixel pitch. A reasonable initial guess can be the centroid + of the cluster. Default is the center of the pixel (0.0, 0.0). + + Returns + ------- + m : Minuit + The minimizer object containing all the information about the fit. + """ # Unpack the grid limits on the x and y axes. xmin, xmax = xlims ymin, ymax = ylims @@ -381,17 +337,19 @@ def nll_grad(x: float, y: float, q: float) -> Tuple[float, float]: # Assign a name to the free parameters. parnames = ["x", "y", "q"] # Initialize the minimizer. + p0 = (*p0, np.sum(pha)) m = Minuit(nll, *p0, grad=nll_grad, name=parnames) # Set the limits for the free parameters. m.limits["x"] = (xmin, xmax) m.limits["y"] = (ymin, ymax) - m.limits["q"] = (np.sum(pha) * 0.5, np.sum(pha) * 1.5) + m.limits["q"] = (np.sum(pha) * 0.5, np.sum(pha) * 1.5) # Set the initial step sizes for the minimizer. We choose half the bin # size as default value. This value is automatically adjusted by the # minimizer during the fit, so it is not critical to choose a very # precise value. m.errors["x"] = bin_size / 2 m.errors["y"] = bin_size / 2 + m.errors["q"] = np.sum(pha) * 0.1 # Define the strategy for the minimizer. We use the higher strategy to # avoid numerical problems with the hessian. For an explanation of the # strategy levels, please refer to the iminuit documentation. @@ -401,4 +359,4 @@ def nll_grad(x: float, y: float, q: float) -> Tuple[float, float]: # Return the minimizer object, which contains the best-fit values and # their uncertainties. It's better to return the whole object to allow # the caller to inspect the fit results and diagnostics. - return m \ No newline at end of file + return m From e2ff53539b6acf2475990c535ae48bfe9338d05f Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Sun, 26 Jul 2026 12:25:41 +0200 Subject: [PATCH 5/5] Minor. --- src/hexsample/clustering.py | 6 +++++- src/hexsample/fileio.py | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py index 8af34f7..b01445f 100644 --- a/src/hexsample/clustering.py +++ b/src/hexsample/clustering.py @@ -52,6 +52,7 @@ class Cluster: adc_to_ev: float pos_recon_algorithm: str recon_pars: Optional[dict] = None + _pulse_height = None def __post_init__(self) -> None: """Small cross check on the dimensions of the arrays passed in the constructor. @@ -67,7 +68,9 @@ def size(self) -> int: def pulse_height(self) -> float: """Return the total pulse height of the cluster. """ - return self.pha.sum() + if self._pulse_height is None: + self._pulse_height = self.pha.sum() + return self._pulse_height def energy(self) -> float: """Return the energy of the cluster in eV. @@ -146,6 +149,7 @@ def mle(self, m = mle(self.pha, equal_noise, position_cal.values, position_cal.bin_size, position_cal.xlims, position_cal.ylims, p0=p0) # Calculate the absolute position of the photon from the fit results. + self._pulse_height = m.values["q"] return self.x[0] + m.values["x"] * pitch, self.y[0] + m.values["y"] * pitch def position(self) -> Tuple[float, float]: diff --git a/src/hexsample/fileio.py b/src/hexsample/fileio.py index dd01322..292cc89 100644 --- a/src/hexsample/fileio.py +++ b/src/hexsample/fileio.py @@ -245,9 +245,9 @@ def _fill_recon_row(row: tables.tableextension.Row, event: ReconEvent) -> None: row["livetime"] = event.livetime #row["roi_size"] = event.roi_size row["cluster_size"] = event.cluster.size() + row["posx"], row["posy"] = event.position() row["adc"] = event.adc() row["energy"] = event.energy() - row["posx"], row["posy"] = event.position() row.append()