Olga Sergijenko / gammalib

Commit 05bbdfcdebf2b1f24879e76d2e33d22bd07cfdbb

Authored by Jürgen Knödlseder
1 parent 35b88c84

Use vector response in unbinned likelihood function (#3203)

Inline Side-by-side
Showing 1 changed file with 17 additions and 9 deletions
src/obs/GObservation.cpp
  View file @05bbdfc
... ... @@ -1341,13 +1341,15 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models,
1341 1341 // Initialise likelihood value
1342 1342 double value = 0.0;
1343 1343  
1344   - // Get number of parameters
1345   - int npars = gradient->size();
  1344 + // Get number of events and parameters
  1345 + int nevents = events()->size();
  1346 + int npars = gradient->size();
1346 1347  
1347 1348 // Allocate some working arrays
1348   - int* inx = new int[npars];
1349   - double* values = new double[npars];
1350   - GVector wrk_grad(npars);
  1349 + int* inx = new int[npars];
  1350 + double* values = new double[npars];
  1351 + GMatrixSparse wrk_matrix(nevents, npars);
  1352 + GVector wrk_grad(npars);
1351 1353  
1352 1354 // Determine Npred value and gradient for this observation
1353 1355 double npred_value = this->npred(models, &wrk_grad);
... ... @@ -1357,20 +1359,26 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models,
1357 1359 *npred += npred_value;
1358 1360 *gradient += wrk_grad;
1359 1361  
  1362 + // Compute model and derivative
  1363 + GVector model_vector = this->model(models, &wrk_matrix);
  1364 +
1360 1365 // Iterate over all events
1361   - for (int i = 0; i < events()->size(); ++i) {
  1366 + for (int i = 0; i < nevents; ++i) {
1362 1367  
1363 1368 // Get event pointer
1364 1369 const GEvent* event = (*events())[i];
1365 1370  
1366   - // Get model and derivative
1367   - double model = this->model(models, *event, &wrk_grad);
  1371 + // Get model value
  1372 + double model = model_vector[i];
1368 1373  
1369 1374 // Skip bin if model is too small (avoids -Inf or NaN gradients)
1370 1375 if (model <= minmod) {
1371 1376 continue;
1372 1377 }
1373 1378  
  1379 + // Extract working gradient multiplied by bin size
  1380 + GVector wrk_grad = wrk_matrix.row(i);
  1381 +
1374 1382 // Create index array of non-zero derivatives and initialise working
1375 1383 // array
1376 1384 int ndev = 0;
... ... @@ -1863,7 +1871,7 @@ double GObservation::npred_kern::eval(const double&amp; x)
1863 1871 * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
1864 1872 * \f]
1865 1873 *
1866   - * where
  1874 + * where
1867 1875 * \f$S(p,E,t)\f$ is the source model,
1868 1876 * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
1869 1877 * \f$p'\f$ is the measured photon direction,
... ...