Justine Devin / gammalib

Commit 79885c6009972b9554221921f7414d7d60b483e2

Authored by ellisowen
1 parent d81007d0

Implements eval_gradients in GModelSpectralGauss - not yet checked

Inline Side-by-side
Showing 2 changed files with 55 additions and 36 deletions
include/GModelSpectralGauss.hpp
  View file @79885c6
... ... @@ -96,7 +96,7 @@ public:
96 96 void mean(const GEnergy& mean);
97 97 GEnergy sigma(void) const;
98 98 void sigma(const GEnergy& sigma);
99   -
  99 + void update_eval_cache(const GEnergy& energy) const;
100 100  
101 101  
102 102 protected:
... ... @@ -104,6 +104,7 @@ protected:
104 104 void init_members(void);
105 105 void copy_members(const GModelSpectralGauss& model);
106 106 void free_members(void);
  107 + void update_eval_cache(const GEnergy& energy) const;
107 108  
108 109 // Energy flux integration kernel
109 110 class eflux_kernel : public GFunction {
... ... @@ -125,6 +126,11 @@ protected:
125 126 GModelPar m_norm; //!< Normalization factor
126 127 GModelPar m_mean; //!< Gaussian mean energy
127 128 GModelPar m_sigma; //!< Gaussian energy width
  129 +
  130 + // Cached members used for pre-computations
  131 + mutable double m_last_norm; //!< Last Normalization value
  132 + mutable double m_last_mean; //!< Last Mean value
  133 + mutable double m_last_sigma; //!< Last Sigma value
128 134 };
129 135  
130 136  
... ...
src/model/GModelSpectralGauss.cpp
  View file @79885c6
... ... @@ -240,8 +240,12 @@ GModelSpectralGauss* GModelSpectralGauss::clone(void) const
240 240 * Evaluates
241 241 *
242 242 * \f[
  243 +<<<<<<< HEAD
243 244 * \frac{dN}{dE} = \frac{\norm}{\sqrt{2\pi}\sigma}
244 245 * \exp(\frac{(E-\bar{E})^2}{2\sigma^2})
  246 +=======
  247 + * \frac{dN}{dE}=\frac{\norm}{\sqrt{2\pi}\sigma}\exp(\frac{-(E-\bar{E})^2}{2\sigma^2})
  248 +>>>>>>> Implements eval_gradients in GModelSpectralGauss - not yet checked
245 249 * \f]
246 250 ***************************************************************************/
247 251 double GModelSpectralGauss::eval(const GEnergy& srcEng,
... ... @@ -272,7 +276,7 @@ double GModelSpectralGauss::eval(const GEnergy&amp; srcEng,
272 276 * This method simply calls the eval() method as no analytical gradients will
273 277 * be computed. See the eval() method for details.
274 278 *
275   - * TODO: Implement this
  279 + * TODO: update docstring.
276 280 ***************************************************************************/
277 281 double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng,
278 282 const GTime& srcTime)
... ... @@ -285,22 +289,25 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
285 289 // Update the evaluation cache
286 290 update_eval_cache(srcEng);
287 291  
288   - // Compute function value
  292 + // Compute function terms
289 293 double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
290   - double term2 = (energy - mean) * (energy - mean) / (2 * sigma * sigma);
291   - double value = term1 * std::exp(- term2);
  294 + double term2 = (1 / sigma) * gammalib::inv_sqrt2pi;
  295 + double term3 = (energy - mean) * (energy - mean) / (2 * sigma * sigma);
  296 + double term4 = (energy - mean) / (sigma * sigma);
  297 + double term5 = (norm / (sigma * sigma)) * gammalib::inv_sqrt2pi;
  298 +
  299 + // Compute function value
  300 + double value = term1 * std::exp(- term3);
  301 +
292 302  
293 303 // Compute partial derivatives with respect to the parameter factor
294   - // values. The partial derivatives with respect to the parameter
295   - // values are obtained by division by the scale factor.
296   - double g_norm = (m_norm.is_free())
297   - ? m_norm.scale() * m_last_power : 0.0;
298   - double g_mean = (m_index.is_free())
299   - ? value * m_index.scale() * std::log(m_last_e_norm) : 0.0;
300   - double g_sigma = (m_ecut.is_free())
301   - ? value * m_last_e_cut / m_ecut.factor_value() : 0.0;
302   -
303   - // Set gradients
  304 + // values (partial differentials were determined analytically).
  305 +
  306 + double g_norm = term2 * std::exp(- term3);
  307 + double g_mean = term1 * term4 * std::exp(- term3);
  308 + double g_sigma = - term5 * std::exp(- term3) * (1 - (2 * term3));
  309 +
  310 + // Set gradients
304 311 m_norm.factor_gradient(g_norm);
305 312 m_mean.factor_gradient(g_mean);
306 313 m_sigma.factor_gradient(g_sigma);
... ... @@ -308,14 +315,14 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
308 315 // Compile option: Check for NaN/Inf
309 316 #if defined(G_NAN_CHECK)
310 317 if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
311   - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients";
  318 + std::cout << "*** ERROR: GModelSpectralGauss::eval_gradients";
312 319 std::cout << "(srcEng=" << srcEng;
313 320 std::cout << ", srcTime=" << srcTime << "):";
314 321 std::cout << " NaN/Inf encountered";
315 322 std::cout << " (value=" << value;
316   - std::cout << ", e_norm=" << m_last_e_norm;
317   - std::cout << ", e_cut=" << m_last_e_cut;
318   - std::cout << ", power=" << m_last_power;
  323 + std::cout << ", norm=" << m_last_norm;
  324 + std::cout << ", mean=" << m_last_mean;
  325 + std::cout << ", Sigma=" << m_last_sigma;
319 326 std::cout << ")" << std::endl;
320 327 }
321 328 #endif
... ... @@ -719,6 +726,11 @@ void GModelSpectralGauss::init_members(void)
719 726 m_sigma.gradient(0.0);
720 727 m_sigma.has_grad(false);
721 728  
  729 + // Initialise eval cache
  730 + m_last_norm = 0.0;
  731 + m_last_mean = 0.0;
  732 + m_last_sigma = 0.0;
  733 +
722 734 // Set parameter pointer(s)
723 735 m_pars.clear();
724 736 m_pars.push_back(&m_norm);
... ... @@ -748,6 +760,11 @@ void GModelSpectralGauss::copy_members(const GModelSpectralGauss&amp; model)
748 760 m_pars.push_back(&m_mean);
749 761 m_pars.push_back(&m_sigma);
750 762  
  763 + // Copy eval cache
  764 + m_last_norm = model.m_last_norm;
  765 + m_last_mean = model.m_last_mean;
  766 + m_last_sigma = model.m_last_sigma;
  767 +
751 768 // Return
752 769 return;
753 770 }
... ... @@ -785,35 +802,31 @@ double GModelSpectralGauss::eflux_kernel::eval(const double&amp; energy)
785 802 *
786 803 * Updates the precomputation cache for eval() and eval_gradients() methods.
787 804 *
788   - * TODO integrate this to correct form
789 805 ***************************************************************************/
790 806 void GModelSpectralGauss::update_eval_cache(const GEnergy& energy) const
791 807 {
792 808 // Get parameter values (takes 3 multiplications which are difficult
793 809 // to avoid)
794   - double index = m_index.value();
795   - double ecut = m_ecut.value();
796   - double pivot = m_pivot.value();
  810 + double norm = m_norm.value();
  811 + double mean = m_mean.value();
  812 + double sigma = m_sigma.value();
797 813  
798   - // If the energy or one of the parameters index, cut-off or pivot
  814 + // If the energy or one of the parameters norm, mean or sigma
799 815 // energy has changed then recompute the cache
800   - if ((m_last_energy != energy) ||
801   - (m_last_index != index) ||
802   - (m_last_ecut != ecut) ||
803   - (m_last_pivot != pivot)) {
  816 + if ((m_last_norm != norm) ||
  817 + (m_last_mean != mean) ||
  818 + (m_last_sigma != sigma)) {
804 819  
805 820 // Store actual energy and parameter values
806   - m_last_energy = energy;
807   - m_last_index = index;
808   - m_last_ecut = ecut;
809   - m_last_pivot = pivot;
  821 + m_last_norm = norm;
  822 + m_last_mean = mean;
  823 + m_last_sigma = sigma;
810 824  
811 825 // Compute and store value
812 826 double eng = energy.MeV();
813   - m_last_e_norm = eng / m_last_pivot;
814   - m_last_e_cut = eng / m_last_ecut;
815   - m_last_power = std::pow(m_last_e_norm, m_last_index) *
816   - std::exp(-m_last_e_cut);
  827 + m_last_norm = eng / m_last_norm;
  828 + m_last_mean = eng / m_last_mean;
  829 + m_last_sigma = eng/ m_last_sigma;
817 830  
818 831 } // endif: recomputation was required
819 832  
... ...