Eirik Hatlen / gammalib

Commit 6b76cc08fc3666043e4b20183d83b69dfbcf7958

Authored by Jürgen Knödlseder
1 parent 34d097cf

Use transposed sparse matrix to store gradients (#3203)

Inline Side-by-side
Showing 3 changed files with 85 additions and 74 deletions
src/model/GModel.cpp
  View file @6b76cc0
... ... @@ -253,17 +253,17 @@ GVector GModel::eval(const GObservation& obs,
253 253  
254 254 // Check matrix consistency
255 255 if (grad) {
256   - if (gradients->rows() != npars) {
257   - std::string msg = "Number of "+gammalib::str(gradients->rows())+
258   - " rows in gradient matrix differs from number "
  256 + if (gradients->columns() != npars) {
  257 + std::string msg = "Number of "+gammalib::str(gradients->columns())+
  258 + " columns in gradient matrix differs from number "
259 259 "of "+gammalib::str(npars)+" parameters "
260 260 "in model. Please provide a compatible gradient "
261 261 "matrix.";
262 262 throw GException::invalid_argument(G_EVAL, msg);
263 263 }
264   - if (gradients->columns() != nevents) {
265   - std::string msg = "Number of "+gammalib::str(gradients->columns())+
266   - " columns in gradient matrix differs from number "
  264 + if (gradients->rows() != nevents) {
  265 + std::string msg = "Number of "+gammalib::str(gradients->rows())+
  266 + " rows in gradient matrix differs from number "
267 267 "of "+gammalib::str(nevents)+" events in "
268 268 "observation. Please provide a compatible "
269 269 "gradient matrix.";
... ... @@ -274,6 +274,12 @@ GVector GModel::eval(const GObservation& obs,
274 274 // Initialise values
275 275 GVector values(nevents);
276 276  
  277 + // Initialise temporary vectors to hold gradients
  278 + GVector* tmp_gradients = new GVector[npars];
  279 + for (int i = 0; i < npars; ++i) {
  280 + tmp_gradients[i] = GVector(nevents);
  281 + }
  282 +
277 283 // Loop over events
278 284 for (int k = 0; k < nevents; ++k) {
279 285  
... ... @@ -287,23 +293,32 @@ GVector GModel::eval(const GObservation&amp; obs,
287 293 if (grad) {
288 294  
289 295 // Initialise gradient vector
290   - GVector gradient(npars);
  296 + //GVector gradient(npars);
291 297  
292 298 // Extract relevant parameter gradients
293 299 for (int i = 0; i < npars; ++i) {
294 300 const GModelPar& par = (*this)[i];
295 301 if (par.is_free() && par.has_grad()) {
296   - gradient[i] = par.factor_gradient();
  302 + //gradient[i] = par.factor_gradient();
  303 + tmp_gradients[i][k] = par.factor_gradient();
297 304 }
298 305 }
299 306  
300 307 // Insert gradient vector into matrix
301   - gradients->add_to_column(k, gradient);
  308 + //gradients->add_to_column(k, gradient);
302 309  
303 310 } // endif: looped over gradients
304 311  
305 312 } // endfor: looped over events
306 313  
  314 + // Fill gradients into matrix
  315 + for (int i = 0; i < npars; ++i) {
  316 + gradients->column(i, tmp_gradients[i]);
  317 + }
  318 +
  319 + // Delete temporal gradients
  320 + delete [] tmp_gradients;
  321 +
307 322 // Return values
308 323 return values;
309 324 }
... ...
src/obs/GObservation.cpp
  View file @6b76cc0
... ... @@ -463,10 +463,12 @@ GVector GObservation::model(const GModels&amp; models,
463 463 else {
464 464  
465 465 // Allocate gradient matrix
466   - GMatrixSparse gradients(mptr->size(), nevents);
  466 + //GMatrixSparse gradients(mptr->size(), nevents);
  467 + GMatrixSparse gradients(nevents, mptr->size());
467 468  
468 469 // Initialise sparse matrix stack
469   - gradients.stack_init(mptr->size()*1000, 1000);
  470 + //gradients.stack_init(mptr->size()*1000, 1000);
  471 + gradients.stack_init(nevents*mptr->size(), mptr->size());
470 472  
471 473 // Evaluate model and add to values
472 474 values += mptr->eval(*this, &gradients);
... ... @@ -495,7 +497,8 @@ GVector GObservation::model(const GModels&amp; models,
495 497 GVector grad(nevents);
496 498 if (par.is_free()) {
497 499 if (par.has_grad()) {
498   - grad = gradients.row(ipar);
  500 + //grad = gradients.row(ipar);
  501 + grad = gradients.column(ipar);
499 502 }
500 503 else {
501 504 grad = model_grad(*mptr, par);
... ... @@ -522,7 +525,8 @@ GVector GObservation::model(const GModels&amp; models,
522 525 GVector grad(nevents);
523 526 if (par.is_free()) {
524 527 if (par.has_grad()) {
525   - grad = gradients.row(ipar);
  528 + //grad = gradients.row(ipar);
  529 + grad = gradients.column(ipar);
526 530 // Debugging
527 531 /*
528 532 if (ipar < 3) {
... ...
src/obs/GResponse.cpp
  View file @6b76cc0
... ... @@ -348,17 +348,17 @@ GVector GResponse::convolve(const GModelSky&amp; model,
348 348  
349 349 // Check matrix consistency
350 350 if (grad) {
351   - if (gradients->rows() != npars) {
352   - std::string msg = "Number of "+gammalib::str(gradients->rows())+
353   - " rows in gradient matrix differs from number "
  351 + if (gradients->columns() != npars) {
  352 + std::string msg = "Number of "+gammalib::str(gradients->columns())+
  353 + " columns in gradient matrix differs from number "
354 354 "of "+gammalib::str(npars)+" parameters "
355 355 "in model. Please provide a compatible gradient "
356 356 "matrix.";
357 357 throw GException::invalid_argument(G_CONVOLVE, msg);
358 358 }
359   - if (gradients->columns() != nevents) {
360   - std::string msg = "Number of "+gammalib::str(gradients->columns())+
361   - " columns in gradient matrix differs from number "
  359 + if (gradients->rows() != nevents) {
  360 + std::string msg = "Number of "+gammalib::str(gradients->rows())+
  361 + " rows in gradient matrix differs from number "
362 362 "of "+gammalib::str(nevents)+" events in "
363 363 "observation. Please provide a compatible "
364 364 "gradient matrix.";
... ... @@ -375,6 +375,12 @@ GVector GResponse::convolve(const GModelSky&amp; model,
375 375 // Case A: Integration
376 376 if (use_edisp()) {
377 377  
  378 + // Initialise temporary vectors to hold gradients
  379 + GVector* tmp_gradients = new GVector[npars];
  380 + for (int i = 0; i < npars; ++i) {
  381 + tmp_gradients[i] = GVector(nevents);
  382 + }
  383 +
378 384 // Loop over events
379 385 for (int k = 0; k < nevents; ++k) {
380 386  
... ... @@ -429,7 +435,7 @@ GVector GResponse::convolve(const GModelSky&amp; model,
429 435 if (grad) {
430 436  
431 437 // Initialise gradient vector
432   - GVector gradient(npars);
  438 + //GVector gradient(npars);
433 439  
434 440 // Get number of spatial, spectral and temporal parameters
435 441 int n_spat = model.spatial()->size();
... ... @@ -442,7 +448,8 @@ GVector GResponse::convolve(const GModelSky&amp; model,
442 448 for (int i = 0; i < n_spec; ++i) {
443 449 GModelPar& par = (*(model.spectral()))[i];
444 450 if (par.is_free() && par.has_grad()) {
445   - gradient[offset+i] = array[index++];
  451 + //gradient[offset+i] = array[index++];
  452 + tmp_gradients[offset+i][k] = array[index++];
446 453 }
447 454 }
448 455 }
... ... @@ -453,7 +460,8 @@ GVector GResponse::convolve(const GModelSky&amp; model,
453 460 for (int i = 0; i < n_temp; ++i) {
454 461 GModelPar& par = (*(model.temporal()))[i];
455 462 if (par.is_free() && par.has_grad()) {
456   - gradient[offset+i] = array[index++];
  463 + //gradient[offset+i] = array[index++];
  464 + tmp_gradients[offset+i][k] = array[index++];
457 465 }
458 466 }
459 467 }
... ... @@ -465,14 +473,15 @@ GVector GResponse::convolve(const GModelSky&amp; model,
465 473 const GModelPar& par = model.scale(i);
466 474 if (par.name() == obs.instrument()) {
467 475 if (par.is_free() && par.has_grad()) {
468   - gradient[offset+i] = array[index++];
  476 + //gradient[offset+i] = array[index++];
  477 + tmp_gradients[offset+i][k] = array[index++];
469 478 }
470 479 }
471 480 }
472 481 }
473 482  
474 483 // Insert gradient vector into matrix
475   - gradients->add_to_column(k, gradient);
  484 + //gradients->add_to_column(k, gradient);
476 485  
477 486 } // endif: set gradients
478 487  
... ... @@ -490,6 +499,14 @@ GVector GResponse::convolve(const GModelSky&amp; model,
490 499  
491 500 } // endfor: looped over events
492 501  
  502 + // Fill gradients into matrix
  503 + for (int i = 0; i < npars; ++i) {
  504 + gradients->column(i, tmp_gradients[i]);
  505 + }
  506 +
  507 + // Delete temporal gradients
  508 + delete [] tmp_gradients;
  509 +
493 510 } // endif: Case A
494 511  
495 512 // Case B: No integration (assume no energy dispersion)
... ... @@ -1624,6 +1641,12 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1624 1641 double scale = (model.has_scales())
1625 1642 ? model.scale(obs.instrument()).value() : 1.0;
1626 1643  
  1644 + // Initialise temporary vectors to hold gradients
  1645 + GVector* tmp_gradients = new GVector[npars];
  1646 + for (int i = 0; i < npars; ++i) {
  1647 + tmp_gradients[i] = GVector(nevents);
  1648 + }
  1649 +
1627 1650 // Loop over events
1628 1651 for (int k = 0; k < nevents; ++k) {
1629 1652  
... ... @@ -1655,7 +1678,7 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1655 1678 if (grad) {
1656 1679  
1657 1680 // Initialise gradient vector
1658   - GVector gradient(npars);
  1681 + //GVector gradient(npars);
1659 1682  
1660 1683 // Get number of spatial, spectral and temporal parameters
1661 1684 int n_spat = model.spatial()->size();
... ... @@ -1666,7 +1689,8 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1666 1689 if (model.spatial() != NULL) {
1667 1690 double fact = spec * temp * scale;
1668 1691 for (int i = 0; i < n_spat; ++i) {
1669   - gradient[i] = spat_gradients(k,i) * fact;
  1692 + //gradient[i] = spat_gradients(k,i) * fact;
  1693 + tmp_gradients[i][k] = spat_gradients(k,i) * fact;
1670 1694 }
1671 1695 }
1672 1696  
... ... @@ -1675,7 +1699,8 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1675 1699 int offset = n_spat;
1676 1700 double fact = spat * temp * scale;
1677 1701 for (int i = 0; i < n_spec; ++i) {
1678   - gradient[offset+i] =
  1702 + //gradient[offset+i] =
  1703 + tmp_gradients[offset+i][k] =
1679 1704 (*(model.spectral()))[i].factor_gradient() * fact;
1680 1705 }
1681 1706 }
... ... @@ -1685,7 +1710,8 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1685 1710 int offset = n_spat + n_spec;
1686 1711 double fact = spat * spec * scale;
1687 1712 for (int i = 0; i < n_temp; ++i) {
1688   - gradient[offset+i] =
  1713 + //gradient[offset+i] =
  1714 + tmp_gradients[offset+i][k] =
1689 1715 (*(model.temporal()))[i].factor_gradient() * fact;
1690 1716 }
1691 1717 }
... ... @@ -1699,61 +1725,19 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1699 1725 double g_scale = 0.0;
1700 1726 if (par.name() == obs.instrument()) {
1701 1727 if (par.is_free()) {
1702   - gradient[offset+i] = par.scale() * fact;
  1728 + //gradient[offset+i] = par.scale() * fact;
  1729 + tmp_gradients[offset+i][k] = par.scale() * fact;
1703 1730 }
1704 1731 }
1705 1732 }
1706 1733 }
1707 1734  
1708 1735 // Add gradient vector to matrix
1709   - gradients->add_to_column(k, gradient);
  1736 + //gradients->add_to_column(k, gradient);
1710 1737  
1711 1738 } // endif: computed partial derivatives
1712 1739  
1713   - } // endif: IRF value was positive
1714   -
1715   - // ... otherwise if gradient computation is requested then set the
1716   - // spectral and temporal gradients to zero
1717   - /*
1718   - else if (grad) {
1719   -
1720   - // Initialise gradient vector
1721   - GVector gradient(npars);
1722   -
1723   - // Get number of spatial, spectral and temporal parameters
1724   - int n_spat = model.spatial()->size();
1725   - int n_spec = model.spectral()->size();
1726   - int n_temp = model.temporal()->size();
1727   -
1728   - // Reset spectral gradients
1729   - if (model.spectral() != NULL) {
1730   - int offset = n_spat;
1731   - for (int i = 0; i < n_spec; ++i) {
1732   - gradient[offset+i] = 0.0;
1733   - }
1734   - }
1735   -
1736   - // Reset temporal gradients
1737   - if (model.temporal() != NULL) {
1738   - int offset = n_spat + n_spec;
1739   - for (int i = 0; i < n_temp; ++i) {
1740   - gradient[offset+i] = 0.0;
1741   - }
1742   - }
1743   -
1744   - // Optionally reset scale gradients
1745   - if (model.has_scales()) {
1746   - int offset = n_spat + n_spec + n_temp;
1747   - for (int i = 0; i < model.scales(); ++i) {
1748   - gradient[offset+i] = 0.0;
1749   - }
1750   - }
1751   -
1752   - // Add gradient vector to matrix
1753   - gradients->add_to_column(k, gradient);
1754   -
1755   - } // endelse: IRF value was not positive
1756   - */
  1740 + } // endif: IRF spatial value was positive
1757 1741  
1758 1742 // Compile option: Check for NaN/Inf
1759 1743 #if defined(G_NAN_CHECK)
... ... @@ -1768,6 +1752,14 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1768 1752  
1769 1753 } // endfor: looped over events
1770 1754  
  1755 + // Fill gradients into matrix
  1756 + for (int i = 0; i < npars; ++i) {
  1757 + gradients->column(i, tmp_gradients[i]);
  1758 + }
  1759 +
  1760 + // Delete temporal gradients
  1761 + delete [] tmp_gradients;
  1762 +
1771 1763 } // endif: Gamma-ray source model had a spatial component
1772 1764  
1773 1765 // Return event probabilities
... ...