Pierrick Martin / gammalib

Commit ed7891e8d0f6365b0a467e4949722ae6a030073a

Authored by Jürgen Knödlseder
1 parent c2540409

Correct parameter gradient computation. 

There were two bugs in the parameter gradient computation.

First, the wrong index was used to compute the break value gradient.

Second, the index gradients were both computed, irrespectively of the energy. However, if the energy is smaller than the break energy, only the first index gradient should be computed while the second will be 0. If the energy is larger, the first index gradient should be 0 and only the second should be computed.

Both bugs have been fixed.
Inline Side-by-side
Showing 2 changed files with 77 additions and 69 deletions
include/GModelSpectralBrokenPlaw.hpp
  View file @ed7891e
... ... @@ -110,27 +110,27 @@ protected:
110 110 // Protected members
111 111 GModelPar m_norm; //!< Normalization factor
112 112 GModelPar m_index1; //!< Spectral index1
  113 + GModelPar m_index2; //!< Spectral index2
113 114 GModelPar m_breakenergy; //!< Energy of spectral break
114   - GModelPar m_index2; //!< Spectral index1
115 115  
116 116 // Cached members used for pre-computations
117   - mutable GEnergy m_last_energy; //!< Last energy value
118   - mutable double m_last_index1; //!< Last index1 parameter
119   - mutable double m_last_index2; //!< Last index1 parameter
120   - mutable double m_last_breakenergy; //!< Last breakenergy parameter
121   - mutable double m_last_e_norm; //!< Last E/Ebreakenergy value
122   - mutable double m_last_log_e_norm; //!< Last ln(E/Ebreakenergy) value
123   - mutable double m_last_power; //!< Last power value
124   - mutable double m_mc_emin; //!< Minimum energy
125   - mutable double m_mc_emax; //!< Maximum energy
126   - mutable double m_mc_exponent1; //!< Exponent (index1+1)
127   - mutable double m_mc_exponent2; //!< Exponent (index2+1)
128   - mutable double m_mc_pow_emin; //!< Power of minimum energy
129   - mutable double m_mc_pow_ewidth; //!< Power of energy width
130   - mutable std::vector<double> m_mc_cum; //!< Cumulative distribution
131   - mutable std::vector<double> m_mc_min; //!< Lower boundary for MC
132   - mutable std::vector<double> m_mc_max; //!< Upper boundary for MC
133   - mutable std::vector<double> m_mc_exp; //!< Exponent for MC
  117 + mutable GEnergy m_last_energy; //!< Last energy value
  118 + mutable double m_last_index1; //!< Last index1 parameter
  119 + mutable double m_last_index2; //!< Last index1 parameter
  120 + mutable double m_last_breakenergy; //!< Last breakenergy parameter
  121 + mutable double m_last_e_norm; //!< Last E/Ebreakenergy value
  122 + mutable double m_last_log_e_norm; //!< Last ln(E/Ebreakenergy) value
  123 + mutable double m_last_power; //!< Last power value
  124 + mutable double m_mc_emin; //!< Minimum energy
  125 + mutable double m_mc_emax; //!< Maximum energy
  126 + mutable double m_mc_exponent1; //!< Exponent (index1+1)
  127 + mutable double m_mc_exponent2; //!< Exponent (index2+1)
  128 + mutable double m_mc_pow_emin; //!< Power of minimum energy
  129 + mutable double m_mc_pow_ewidth; //!< Power of energy width
  130 + mutable std::vector<double> m_mc_cum; //!< Cumulative distribution
  131 + mutable std::vector<double> m_mc_min; //!< Lower boundary for MC
  132 + mutable std::vector<double> m_mc_max; //!< Upper boundary for MC
  133 + mutable std::vector<double> m_mc_exp; //!< Exponent for MC
134 134 };
135 135  
136 136  
... ...
src/model/GModelSpectralBrokenPlaw.cpp
  View file @ed7891e
... ... @@ -338,24 +338,34 @@ double GModelSpectralBrokenPlaw::eval_gradients(const GEnergy&amp; srcEng,
338 338 // Compute function value
339 339 double value = m_norm.value() * m_last_power;
340 340  
341   - // Compute partial derivatives of the parameter values
  341 + // Compute normalisation gradient
342 342 double g_norm = (m_norm.isfree())
343 343 ? m_norm.scale() * m_last_power : 0.0;
344   - double g_index1 = (m_index1.isfree())
345   - ? value * m_index1.scale() * m_last_log_e_norm : 0.0;
346   - double g_index2 = (m_index2.isfree())
347   - ? value * m_index2.scale() * m_last_log_e_norm : 0.0;
348   - double index_current = (srcEng.MeV() > m_breakenergy.value())
349   - ? m_last_index1 : m_last_index2;
350   - double g_breakenergy = (m_breakenergy.isfree())
351   - ? -value * index_current / m_breakenergy.factor_value() : 0.0;
352   -
353   -
354   - // Set gradients
355 344 m_norm.factor_gradient(g_norm);
356   - m_index1.factor_gradient(g_index1);
357   - m_breakenergy.factor_gradient(g_breakenergy);
358   - m_index2.factor_gradient(g_index2);
  345 +
  346 + // Compute index and break value gradients
  347 + if (srcEng.MeV() < m_breakenergy.value()) {
  348 + double g_index = (m_index1.isfree())
  349 + ? value * m_index1.scale() * m_last_log_e_norm
  350 + : 0.0;
  351 + double g_break = (m_breakenergy.isfree())
  352 + ? -value * m_last_index1 / m_breakenergy.factor_value()
  353 + : 0.0;
  354 + m_index1.factor_gradient(g_index);
  355 + m_index2.factor_gradient(0.0);
  356 + m_breakenergy.factor_gradient(g_break);
  357 + }
  358 + else {
  359 + double g_index = (m_index2.isfree())
  360 + ? value * m_index2.scale() * m_last_log_e_norm
  361 + : 0.0;
  362 + double g_break = (m_breakenergy.isfree())
  363 + ? -value * m_last_index2 / m_breakenergy.factor_value()
  364 + : 0.0;
  365 + m_index1.factor_gradient(0.0);
  366 + m_index2.factor_gradient(g_index);
  367 + m_breakenergy.factor_gradient(g_break);
  368 + }
359 369  
360 370 // Compile option: Check for NaN/Inf
361 371 #if defined(G_NAN_CHECK)
... ... @@ -368,9 +378,6 @@ double GModelSpectralBrokenPlaw::eval_gradients(const GEnergy&amp; srcEng,
368 378 std::cout << ", e_norm=" << m_last_e_norm;
369 379 std::cout << ", power=" << m_last_power;
370 380 std::cout << ", g_norm=" << g_norm;
371   - std::cout << ", g_index1=" << g_index1;
372   - std::cout << ", g_breakenergy=" << g_breakenergy;
373   - std::cout << ", g_index2=" << g_index2;
374 381 std::cout << ")" << std::endl;
375 382 }
376 383 #endif
... ... @@ -806,16 +813,6 @@ void GModelSpectralBrokenPlaw::init_members(void)
806 813 m_index1.gradient(0.0);
807 814 m_index1.hasgrad(true);
808 815  
809   - // Initialise break energy
810   - m_breakenergy.clear();
811   - m_breakenergy.name("BreakValue");
812   - m_breakenergy.unit("MeV");
813   - m_breakenergy.scale(1.0);
814   - m_breakenergy.value(100.0); // default: 100
815   - m_breakenergy.fix();
816   - m_breakenergy.gradient(0.0);
817   - m_breakenergy.hasgrad(true);
818   -
819 816 // Initialise powerlaw index2
820 817 m_index2.clear();
821 818 m_index2.name("Index2");
... ... @@ -826,6 +823,16 @@ void GModelSpectralBrokenPlaw::init_members(void)
826 823 m_index2.gradient(0.0);
827 824 m_index2.hasgrad(true);
828 825  
  826 + // Initialise break energy
  827 + m_breakenergy.clear();
  828 + m_breakenergy.name("BreakValue");
  829 + m_breakenergy.unit("MeV");
  830 + m_breakenergy.scale(1.0);
  831 + m_breakenergy.value(100.0); // default: 100
  832 + m_breakenergy.fix();
  833 + m_breakenergy.gradient(0.0);
  834 + m_breakenergy.hasgrad(true);
  835 +
829 836 // Set parameter pointer(s)
830 837 m_pars.clear();
831 838 m_pars.push_back(&m_norm);
... ... @@ -925,8 +932,8 @@ void GModelSpectralBrokenPlaw::update_eval_cache(const GEnergy&amp; energy) const
925 932 double breakenergy = m_breakenergy.value();
926 933 double index2 = m_index2.value();
927 934  
928   - // If the energy or one of the parameters index1, index2 or breakenergy energy has
929   - // changed then recompute the cache
  935 + // If the energy or one of the parameters index1, index2 or breakenergy
  936 + // energy has changed then recompute the cache
930 937 if ((m_last_energy != energy) ||
931 938 (m_last_index1 != index1) ||
932 939 (m_last_index2 != index2) ||
... ... @@ -942,7 +949,6 @@ void GModelSpectralBrokenPlaw::update_eval_cache(const GEnergy&amp; energy) const
942 949 double eng = energy.MeV();
943 950 m_last_e_norm = eng / m_last_breakenergy;
944 951 m_last_log_e_norm = std::log(m_last_e_norm);
945   -
946 952 if (eng < m_last_breakenergy ) {
947 953 m_last_power = std::pow(m_last_e_norm, m_last_index1);
948 954 }
... ... @@ -966,9 +972,9 @@ void GModelSpectralBrokenPlaw::update_eval_cache(const GEnergy&amp; energy) const
966 972 * Updates the precomputation cache for Monte Carlo simulations.
967 973 ***************************************************************************/
968 974 void GModelSpectralBrokenPlaw::update_mc_cache(const GEnergy& emin,
969   - const GEnergy& emax) const
970   -
971   -{ // Check if we need to update the cache
  975 + const GEnergy& emax) const
  976 +{
  977 + // Check if we need to update the cache
972 978 if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
973 979  
974 980 // Store new energy interval
... ... @@ -981,7 +987,6 @@ void GModelSpectralBrokenPlaw::update_mc_cache(const GEnergy&amp; emin,
981 987 m_mc_max.clear();
982 988 m_mc_exp.clear();
983 989  
984   -
985 990 // Get energy range in MeV
986 991 double e_min = emin.MeV();
987 992 double e_max = emax.MeV();
... ... @@ -992,19 +997,20 @@ void GModelSpectralBrokenPlaw::update_mc_cache(const GEnergy&amp; emin,
992 997 // Allocate flux
993 998 double flux;
994 999  
995   - // Determine left node index for minimum energy
  1000 + // Determine left node index for minimum and maximum energy
996 1001 int inx_emin = (e_min < m_breakenergy.value() ) ? 0 : 1;
997   - // Determine left node index for maximum energy
998 1002 int inx_emax = (e_max < m_breakenergy.value() ) ? 0 : 1;
999 1003  
1000 1004 // If both energies are within the same node then just
1001 1005 // add this one node on the stack
1002 1006 if (inx_emin == inx_emax) {
1003   - double exp_valid = (e_min < m_breakenergy.value()) ? m_index1.value() : m_index2.value();
1004   - flux = m_norm.value()*gammalib::plaw_photon_flux(e_min,
1005   - e_max,
1006   - m_breakenergy.value(),
1007   - exp_valid);
  1007 + double exp_valid = (e_min < m_breakenergy.value())
  1008 + ? m_index1.value() : m_index2.value();
  1009 + flux = m_norm.value() *
  1010 + gammalib::plaw_photon_flux(e_min,
  1011 + e_max,
  1012 + m_breakenergy.value(),
  1013 + exp_valid);
1008 1014 m_mc_cum.push_back(flux);
1009 1015 m_mc_min.push_back(e_min);
1010 1016 m_mc_max.push_back(e_max);
... ... @@ -1014,27 +1020,28 @@ void GModelSpectralBrokenPlaw::update_mc_cache(const GEnergy&amp; emin,
1014 1020 // ... otherwise integrate over both nodes
1015 1021 else {
1016 1022 // just enter the values for first pl: bin [0]
1017   - flux = m_norm.value()*gammalib::plaw_photon_flux(e_min,
1018   - m_breakenergy.value(),
1019   - m_breakenergy.value(),
1020   - m_index1.value());
  1023 + flux = m_norm.value() *
  1024 + gammalib::plaw_photon_flux(e_min,
  1025 + m_breakenergy.value(),
  1026 + m_breakenergy.value(),
  1027 + m_index1.value());
1021 1028 m_mc_cum.push_back(flux);
1022 1029 m_mc_exp.push_back(m_index1.value());
1023 1030 m_mc_min.push_back(e_min);
1024 1031 m_mc_max.push_back(m_breakenergy.value());
1025 1032  
1026 1033 // and for
1027   - flux = m_norm.value()*gammalib::plaw_photon_flux(m_breakenergy.value(),
1028   - e_max,
1029   - m_breakenergy.value(),
1030   - m_index2.value());
  1034 + flux = m_norm.value() *
  1035 + gammalib::plaw_photon_flux(m_breakenergy.value(),
  1036 + e_max,
  1037 + m_breakenergy.value(),
  1038 + m_index2.value());
1031 1039 m_mc_cum.push_back(flux);
1032 1040 m_mc_exp.push_back(m_index2.value());
1033 1041 m_mc_max.push_back(e_max);
1034 1042 m_mc_min.push_back(m_breakenergy.value());
1035 1043 } // endelse: emin and emax not between same nodes
1036 1044  
1037   -
1038 1045 // Build cumulative distribution
1039 1046 for (int i = 1; i < m_mc_cum.size(); ++i) {
1040 1047 m_mc_cum[i] += m_mc_cum[i-1];
... ... @@ -1046,6 +1053,7 @@ void GModelSpectralBrokenPlaw::update_mc_cache(const GEnergy&amp; emin,
1046 1053  
1047 1054 // Set MC values
1048 1055 for (int i = 0; i < m_mc_cum.size(); ++i) {
  1056 +
1049 1057 // Compute exponent
1050 1058 double exponent = m_mc_exp[i] + 1.0;
1051 1059  
... ...