From f0265fbbe9a6fb74e45fda96aad3df4e6c20ec5a Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 15 Jan 2021 15:52:18 +0100 Subject: [PATCH] Transpose matrix in GObservation::likelihood methods for faster column-wise instead of row-wise access (#3501) --- NEWS | 4 ++++ src/obs/GObservation.cpp | 28 ++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+), 0 deletions(-) diff --git a/NEWS b/NEWS index e1ce116..1cf70b7 100644 --- a/NEWS +++ b/NEWS @@ -136,6 +136,10 @@ of functions (#3328). 10. Observation module ---------------------- +The sparse gradient matrix is now transposed in the likelihood method so that the +elements can be access column-wise instead of row-wise, which speeds-up the element +access (#3501). + Add source removal from vector cache to GResponse::remove_response_cache() method (#3499). diff --git a/src/obs/GObservation.cpp b/src/obs/GObservation.cpp index afd51d4..d4ca212 100644 --- a/src/obs/GObservation.cpp +++ b/src/obs/GObservation.cpp @@ -62,6 +62,7 @@ const double minerr = 1.0e-100; //!< Minimum statistical error /* __ Coding definitions _________________________________________________ */ #define G_LN_ENERGY_INT //!< ln(E) variable substitution for integration +#define G_TRANSPOSE_MATRIX //!< Transpose matrix in likelihood methods /* __ Debug definitions __________________________________________________ */ //#define G_OPT_DEBUG //!< Debug likelihood computation @@ -1394,6 +1395,11 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models, // Compute model and derivative GVector model_vector = this->model(models, &wrk_matrix); + // Compile option: transpose matrix + #if defined(G_TRANSPOSE_MATRIX) + wrk_matrix = wrk_matrix.transpose(); + #endif + // Iterate over all events for (int i = 0; i < nevents; ++i) { @@ -1409,7 +1415,11 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models, } // Extract working gradient multiplied by bin size + #if defined(G_TRANSPOSE_MATRIX) + GVector wrk_grad = wrk_matrix.column(i); + #else GVector wrk_grad = wrk_matrix.row(i); + #endif // Create index array of non-zero derivatives and initialise working // array @@ -1525,6 +1535,11 @@ double GObservation::likelihood_poisson_binned(const GModels& models, // Compute model and derivative GVector model_vector = this->model(models, &wrk_matrix); + // Compile option: transpose matrix + #if defined(G_TRANSPOSE_MATRIX) + wrk_matrix = wrk_matrix.transpose(); + #endif + // Iterate over all bins for (int i = 0; i < nevents; ++i) { @@ -1570,7 +1585,11 @@ double GObservation::likelihood_poisson_binned(const GModels& models, *npred += model; // Extract working gradient multiplied by bin size + #if defined(G_TRANSPOSE_MATRIX) + GVector wrk_grad = wrk_matrix.column(i) * bin->size(); + #else GVector wrk_grad = wrk_matrix.row(i) * bin->size(); + #endif // Create index array of non-zero derivatives and initialise working // array @@ -1724,6 +1743,11 @@ double GObservation::likelihood_gaussian_binned(const GModels& models, // Compute model and derivative GVector model_vector = this->model(models, &wrk_matrix); + // Compile option: transpose matrix + #if defined(G_TRANSPOSE_MATRIX) + wrk_matrix = wrk_matrix.transpose(); + #endif + // Iterate over all bins for (int i = 0; i < nevents; ++i) { @@ -1759,7 +1783,11 @@ double GObservation::likelihood_gaussian_binned(const GModels& models, *npred += model; // Extract working gradient multiplied by bin size + #if defined(G_TRANSPOSE_MATRIX) + GVector wrk_grad = wrk_matrix.column(i) * bin->size(); + #else GVector wrk_grad = wrk_matrix.row(i) * bin->size(); + #endif // Create index array of non-zero derivatives and initialise working // array -- libgit2 0.21.4