Eirik Hatlen / gammalib

Commit f0265fbbe9a6fb74e45fda96aad3df4e6c20ec5a

Authored by Jürgen Knödlseder
1 parent 00269f87

Transpose matrix in GObservation::likelihood methods for faster column-wise inst… 

…ead of row-wise access (#3501)
Inline Side-by-side
Showing 2 changed files with 32 additions and 0 deletions
... ... @@ -136,6 +136,10 @@ of functions (#3328).
136 136  
137 137 10. Observation module
138 138 ----------------------
  139 +The sparse gradient matrix is now transposed in the likelihood method so that the
  140 +elements can be access column-wise instead of row-wise, which speeds-up the element
  141 +access (#3501).
  142 +
139 143 Add source removal from vector cache to GResponse::remove_response_cache()
140 144 method (#3499).
141 145  
... ...
src/obs/GObservation.cpp
  View file @f0265fb
... ... @@ -62,6 +62,7 @@ const double minerr = 1.0e-100; //!< Minimum statistical error
62 62  
63 63 /* __ Coding definitions _________________________________________________ */
64 64 #define G_LN_ENERGY_INT //!< ln(E) variable substitution for integration
  65 +#define G_TRANSPOSE_MATRIX //!< Transpose matrix in likelihood methods
65 66  
66 67 /* __ Debug definitions __________________________________________________ */
67 68 //#define G_OPT_DEBUG //!< Debug likelihood computation
... ... @@ -1394,6 +1395,11 @@ double GObservation::likelihood_poisson_unbinned(const GModels&amp; models,
1394 1395 // Compute model and derivative
1395 1396 GVector model_vector = this->model(models, &wrk_matrix);
1396 1397  
  1398 + // Compile option: transpose matrix
  1399 + #if defined(G_TRANSPOSE_MATRIX)
  1400 + wrk_matrix = wrk_matrix.transpose();
  1401 + #endif
  1402 +
1397 1403 // Iterate over all events
1398 1404 for (int i = 0; i < nevents; ++i) {
1399 1405  
... ... @@ -1409,7 +1415,11 @@ double GObservation::likelihood_poisson_unbinned(const GModels&amp; models,
1409 1415 }
1410 1416  
1411 1417 // Extract working gradient multiplied by bin size
  1418 + #if defined(G_TRANSPOSE_MATRIX)
  1419 + GVector wrk_grad = wrk_matrix.column(i);
  1420 + #else
1412 1421 GVector wrk_grad = wrk_matrix.row(i);
  1422 + #endif
1413 1423  
1414 1424 // Create index array of non-zero derivatives and initialise working
1415 1425 // array
... ... @@ -1525,6 +1535,11 @@ double GObservation::likelihood_poisson_binned(const GModels&amp; models,
1525 1535 // Compute model and derivative
1526 1536 GVector model_vector = this->model(models, &wrk_matrix);
1527 1537  
  1538 + // Compile option: transpose matrix
  1539 + #if defined(G_TRANSPOSE_MATRIX)
  1540 + wrk_matrix = wrk_matrix.transpose();
  1541 + #endif
  1542 +
1528 1543 // Iterate over all bins
1529 1544 for (int i = 0; i < nevents; ++i) {
1530 1545  
... ... @@ -1570,7 +1585,11 @@ double GObservation::likelihood_poisson_binned(const GModels&amp; models,
1570 1585 *npred += model;
1571 1586  
1572 1587 // Extract working gradient multiplied by bin size
  1588 + #if defined(G_TRANSPOSE_MATRIX)
  1589 + GVector wrk_grad = wrk_matrix.column(i) * bin->size();
  1590 + #else
1573 1591 GVector wrk_grad = wrk_matrix.row(i) * bin->size();
  1592 + #endif
1574 1593  
1575 1594 // Create index array of non-zero derivatives and initialise working
1576 1595 // array
... ... @@ -1724,6 +1743,11 @@ double GObservation::likelihood_gaussian_binned(const GModels&amp; models,
1724 1743 // Compute model and derivative
1725 1744 GVector model_vector = this->model(models, &wrk_matrix);
1726 1745  
  1746 + // Compile option: transpose matrix
  1747 + #if defined(G_TRANSPOSE_MATRIX)
  1748 + wrk_matrix = wrk_matrix.transpose();
  1749 + #endif
  1750 +
1727 1751 // Iterate over all bins
1728 1752 for (int i = 0; i < nevents; ++i) {
1729 1753  
... ... @@ -1759,7 +1783,11 @@ double GObservation::likelihood_gaussian_binned(const GModels&amp; models,
1759 1783 *npred += model;
1760 1784  
1761 1785 // Extract working gradient multiplied by bin size
  1786 + #if defined(G_TRANSPOSE_MATRIX)
  1787 + GVector wrk_grad = wrk_matrix.column(i) * bin->size();
  1788 + #else
1762 1789 GVector wrk_grad = wrk_matrix.row(i) * bin->size();
  1790 + #endif
1763 1791  
1764 1792 // Create index array of non-zero derivatives and initialise working
1765 1793 // array
... ...