Jürgen Knödlseder / gammalib

Commit 8f77d64a790ee634e98ce02bc4d9dc2d88aa71cd

Authored by Jürgen Knödlseder
1 parent 24dbc9ae

Compute spectral model gradients for negative factor values

Inline Side-by-side
Showing 13 changed files with 24 additions and 18 deletions
ChangeLog
  View file @8f77d64
1   -2021-11-28
  1 +2021-11-30
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
5 5  
  6 + Compute spectral model gradients for negative factor values
6 7 Add likelihood fit event selection hook (#3912)
7 8 Add GCOMModelDRM class (#3910)
8 9 Add angle(GVector&,GVector&) friend function (#3908)
... ...
1 1 New Features and Important Changes in GammaLib 2.0.0
2 2  
3   -28 November 2021
  3 +30 November 2021
4 4  
5 5  
6 6 1. Introduction
... ... @@ -197,6 +197,10 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods.
197 197  
198 198 9. Model module
199 199 ---------------
  200 +Gradients for spectral model parameters are now also computed for all models in
  201 +case that the factor values are negative. This provides now error bars for photon
  202 +and energy flux power laws that are allowed to take negative flux values.
  203 +
200 204 Added GModelSpectralFunc spectral model constructor that enables creating a file
201 205 function from an arbitrary spectral model (#3626).
202 206  
... ...
doc/source/admin/release_history/2.0.rst
  View file @8f77d64
1 1 .. _2.0:
2 2  
3   -GammaLib 2.0 release (28 November 2021)
4   -==================================
  3 +GammaLib 2.0 release (30 November 2021)
  4 +=======================================
5 5  
6 6 Overview
7 7 --------
... ... @@ -18,6 +18,7 @@ In particular, this release provides:
18 18 Bug fixes
19 19 ---------
20 20  
  21 +* Compute spectral model gradients for negative factor values
21 22 * [`3563 <https://cta-redmine.irap.omp.eu/issues/3563>`_] -
22 23 Fix incompatibility in ``GMatrixSparse::mix_column_prepare()`` method
23 24 * [`3553 <https://cta-redmine.irap.omp.eu/issues/3553>`_] -
... ...
src/model/GModelSpectralBrokenPlaw.cpp
  View file @8f77d64
... ... @@ -367,7 +367,7 @@ double GModelSpectralBrokenPlaw::eval(const GEnergy&amp; srcEng,
367 367 double g_index = (m_index1.is_free())
368 368 ? value * m_index1.scale() * m_last_log_e_norm
369 369 : 0.0;
370   - double g_break = (m_breakenergy.is_free())
  370 + double g_break = (m_breakenergy.is_free() && m_breakenergy.factor_value() != 0.0)
371 371 ? -value * m_last_index1 / m_breakenergy.factor_value()
372 372 : 0.0;
373 373 m_index1.factor_gradient(g_index);
... ... @@ -378,7 +378,7 @@ double GModelSpectralBrokenPlaw::eval(const GEnergy&amp; srcEng,
378 378 double g_index = (m_index2.is_free())
379 379 ? value * m_index2.scale() * m_last_log_e_norm
380 380 : 0.0;
381   - double g_break = (m_breakenergy.is_free())
  381 + double g_break = (m_breakenergy.is_free() && m_breakenergy.factor_value() != 0.0)
382 382 ? -value * m_last_index2 / m_breakenergy.factor_value()
383 383 : 0.0;
384 384 m_index1.factor_gradient(0.0);
... ...
src/model/GModelSpectralExpInvPlaw.cpp
  View file @8f77d64
... ... @@ -388,7 +388,7 @@ double GModelSpectralExpInvPlaw::eval(const GEnergy&amp; srcEng,
388 388 double g_lambda = (m_lambda.is_free())
389 389 ? -value * m_lambda.scale() * m_last_energy.MeV()
390 390 : 0.0;
391   - double g_pivot = (m_pivot.is_free())
  391 + double g_pivot = (m_pivot.is_free() && m_pivot.factor_value() != 0.0)
392 392 ? -value * m_last_index / m_pivot.factor_value()
393 393 : 0.0;
394 394  
... ...
src/model/GModelSpectralExpPlaw.cpp
  View file @8f77d64
... ... @@ -355,10 +355,10 @@ double GModelSpectralExpPlaw::eval(const GEnergy&amp; srcEng,
355 355 double g_index = (m_index.is_free())
356 356 ? value * m_index.scale() * std::log(m_last_e_norm)
357 357 : 0.0;
358   - double g_ecut = (m_ecut.is_free())
  358 + double g_ecut = (m_ecut.is_free() && m_ecut.factor_value() != 0.0)
359 359 ? value * m_last_e_cut / m_ecut.factor_value()
360 360 : 0.0;
361   - double g_pivot = (m_pivot.is_free())
  361 + double g_pivot = (m_pivot.is_free() && m_pivot.factor_value() != 0.0)
362 362 ? -value * m_last_index / m_pivot.factor_value()
363 363 : 0.0;
364 364  
... ...
src/model/GModelSpectralLogParabola.cpp
  View file @8f77d64
... ... @@ -361,7 +361,7 @@ double GModelSpectralLogParabola::eval(const GEnergy&amp; srcEng,
361 361 double g_curvature = (m_curvature.is_free())
362 362 ? value * m_curvature.scale() * m_last_log_e_norm *
363 363 m_last_log_e_norm : 0.0;
364   - double g_pivot = (m_pivot.is_free())
  364 + double g_pivot = (m_pivot.is_free() && m_pivot.factor_value() != 0.0)
365 365 ? -value * (m_last_exponent + m_curvature.value() *
366 366 m_last_log_e_norm) / m_pivot.factor_value() : 0.0;
367 367  
... ...
src/model/GModelSpectralNodes.cpp
  View file @8f77d64
... ... @@ -308,13 +308,13 @@ double GModelSpectralNodes::eval(const GEnergy&amp; srcEng,
308 308 }
309 309  
310 310 // Gradient for left node
311   - if (m_values[inx_left].is_free()) {
  311 + if (m_values[inx_left].is_free() && m_values[inx_left].factor_value() != 0.0) {
312 312 double grad = value * wgt_left / m_values[inx_left].factor_value();
313 313 m_values[inx_left].factor_gradient(grad);
314 314 }
315 315  
316 316 // Gradient for right node
317   - if (m_values[inx_right].is_free()) {
  317 + if (m_values[inx_right].is_free() && m_values[inx_right].factor_value() != 0.0) {
318 318 double grad = value * wgt_right / m_values[inx_right].factor_value();
319 319 m_values[inx_right].factor_gradient(grad);
320 320 }
... ...
src/model/GModelSpectralPlaw.cpp
  View file @8f77d64
... ... @@ -331,7 +331,7 @@ double GModelSpectralPlaw::eval(const GEnergy&amp; srcEng,
331 331 ? m_norm.scale() * m_last_power : 0.0;
332 332 double g_index = (m_index.is_free())
333 333 ? value * m_index.scale() * m_last_log_e_norm : 0.0;
334   - double g_pivot = (m_pivot.is_free())
  334 + double g_pivot = (m_pivot.is_free() && m_pivot.factor_value() != 0.0)
335 335 ? -value * m_last_index / m_pivot.factor_value() : 0.0;
336 336  
337 337 // Set gradients
... ...
src/model/GModelSpectralPlawEnergyFlux.cpp
  View file @8f77d64
... ... @@ -352,7 +352,7 @@ double GModelSpectralPlawEnergyFlux::eval(const GEnergy&amp; srcEng,
352 352 if (gradients) {
353 353  
354 354 // Compute partial derivatives of the parameter values
355   - double g_eflux = (m_eflux.is_free() && m_eflux.factor_value() > 0.0)
  355 + double g_eflux = (m_eflux.is_free() && m_eflux.factor_value() != 0.0)
356 356 ? g_eflux = value / m_eflux.factor_value()
357 357 : 0.0;
358 358 double g_index = (m_index.is_free())
... ...
src/model/GModelSpectralPlawPhotonFlux.cpp
  View file @8f77d64
... ... @@ -360,7 +360,7 @@ double GModelSpectralPlawPhotonFlux::eval(const GEnergy&amp; srcEng,
360 360 if (gradients) {
361 361  
362 362 // Compute partial derivatives of the parameter values
363   - double g_flux = (m_flux.is_free() && m_flux.factor_value() > 0.0)
  363 + double g_flux = (m_flux.is_free() && m_flux.factor_value() != 0.0)
364 364 ? value / m_flux.factor_value()
365 365 : 0.0;
366 366 double g_index = (m_index.is_free())
... ...
src/model/GModelSpectralSmoothBrokenPlaw.cpp
  View file @8f77d64
... ... @@ -409,7 +409,7 @@ double GModelSpectralSmoothBrokenPlaw::eval(const GEnergy&amp; srcEng,
409 409 : 0.0;
410 410  
411 411 // Compute pivot and break energy value gradients
412   - double g_pivot = (m_pivot.is_free())
  412 + double g_pivot = (m_pivot.is_free() && m_pivot.factor_value() != 0.0)
413 413 ? -value * m_last_index1 / m_pivot.factor_value()
414 414 : 0.0;
415 415 double g_break = (m_breakenergy.is_free())
... ...
src/model/GModelSpectralSuperExpPlaw.cpp
  View file @8f77d64
... ... @@ -374,11 +374,11 @@ double GModelSpectralSuperExpPlaw::eval(const GEnergy&amp; srcEng,
374 374 double g_index1 = (m_index1.is_free())
375 375 ? value * m_index1.scale() * std::log(m_last_e_norm)
376 376 : 0.0;
377   - double g_ecut = (m_ecut.is_free())
  377 + double g_ecut = (m_ecut.is_free() && m_ecut.factor_value() != 0.0)
378 378 ? value * m_last_index2 * m_last_exponent /
379 379 m_ecut.factor_value()
380 380 : 0.0;
381   - double g_pivot = (m_pivot.is_free())
  381 + double g_pivot = (m_pivot.is_free() && m_pivot.factor_value() != 0.0)
382 382 ? -value * m_last_index1 / m_pivot.factor_value()
383 383 : 0.0;
384 384 double g_index2 = (m_index2.is_free())
... ...