Johan Bregeon / gammalib

Commit b84e41a8e764ed6aa1324c00aeb8b7639fe23345

Authored by ellisowen
1 parent 94cf8055

Implements GModelSpectralGauss eval_gradients method and changes the GModelSpect… 

…ralGauss mc method as suggested
Inline Side-by-side
Showing 2 changed files with 4 additions and 73 deletions
include/GModelSpectralGauss.hpp
  View file @b84e41a
... ... @@ -127,7 +127,6 @@ protected:
127 127 GModelPar m_mean; //!< Gaussian mean energy
128 128 GModelPar m_sigma; //!< Gaussian energy width
129 129  
130   - // Cached members used for pre-computations
131 130 mutable double m_last_norm; //!< Last Normalization value
132 131 mutable double m_last_mean; //!< Last Mean value
133 132 mutable double m_last_sigma; //!< Last Sigma value
... ... @@ -189,7 +188,7 @@ GEnergy GModelSpectralGauss::mean(void) const
189 188 /***********************************************************************//**
190 189 * @brief Set mean energy
191 190 *
192   - * @param[in] mean Mean energy.
  191 + * @param[in] Mean energy.
193 192 ***************************************************************************/
194 193 inline
195 194 void GModelSpectralGauss::mean(const GEnergy& mean)
... ... @@ -215,7 +214,7 @@ GEnergy GModelSpectralGauss::sigma(void) const
215 214 /***********************************************************************//**
216 215 * @brief Set energy width
217 216 *
218   - * @param[in] sigma Energy width
  217 + * @param[in] Energy width
219 218 ***************************************************************************/
220 219 inline
221 220 void GModelSpectralGauss::sigma(const GEnergy& sigma)
... ... @@ -224,4 +223,5 @@ void GModelSpectralGauss::sigma(const GEnergy&amp; sigma)
224 223 return;
225 224 }
226 225  
  226 +
227 227 #endif /* GMODELSPECTRALGAUSS_HPP */
... ...
src/model/GModelSpectralGauss.cpp
  View file @b84e41a
... ... @@ -240,12 +240,7 @@ GModelSpectralGauss* GModelSpectralGauss::clone(void) const
240 240 * Evaluates
241 241 *
242 242 * \f[
243   -<<<<<<< HEAD
244   - * \frac{dN}{dE} = \frac{\norm}{\sqrt{2\pi}\sigma}
245   - * \exp(\frac{(E-\bar{E})^2}{2\sigma^2})
246   -=======
247 243 * \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
249 244 * \f]
250 245 ***************************************************************************/
251 246 double GModelSpectralGauss::eval(const GEnergy& srcEng,
... ... @@ -276,11 +271,7 @@ double GModelSpectralGauss::eval(const GEnergy&amp; srcEng,
276 271 * This method simply calls the eval() method as no analytical gradients will
277 272 * be computed. See the eval() method for details.
278 273 *
279   -<<<<<<< HEAD
280 274 * TODO: update docstring.
281   -=======
282   - * TODO: Implement this
283   ->>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
284 275 ***************************************************************************/
285 276 double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng,
286 277 const GTime& srcTime)
... ... @@ -293,7 +284,6 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
293 284 // Update the evaluation cache
294 285 update_eval_cache(srcEng);
295 286  
296   -<<<<<<< HEAD
297 287 // Compute function terms
298 288 double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
299 289 double term2 = (1 / sigma) * gammalib::inv_sqrt2pi;
... ... @@ -313,24 +303,6 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
313 303 double g_sigma = - term5 * std::exp(- term3) * (1 - (2 * term3));
314 304  
315 305 // Set gradients
316   -=======
317   - // Compute function value
318   - double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
319   - double term2 = (energy - mean) * (energy - mean) / (2 * sigma * sigma);
320   - double value = term1 * std::exp(- term2);
321   -
322   - // Compute partial derivatives with respect to the parameter factor
323   - // values. The partial derivatives with respect to the parameter
324   - // values are obtained by division by the scale factor.
325   - double g_norm = (m_norm.is_free())
326   - ? m_norm.scale() * m_last_power : 0.0;
327   - double g_mean = (m_index.is_free())
328   - ? value * m_index.scale() * std::log(m_last_e_norm) : 0.0;
329   - double g_sigma = (m_ecut.is_free())
330   - ? value * m_last_e_cut / m_ecut.factor_value() : 0.0;
331   -
332   - // Set gradients
333   ->>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
334 306 m_norm.factor_gradient(g_norm);
335 307 m_mean.factor_gradient(g_mean);
336 308 m_sigma.factor_gradient(g_sigma);
... ... @@ -338,24 +310,14 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
338 310 // Compile option: Check for NaN/Inf
339 311 #if defined(G_NAN_CHECK)
340 312 if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
341   -<<<<<<< HEAD
342 313 std::cout << "*** ERROR: GModelSpectralGauss::eval_gradients";
343   -=======
344   - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients";
345   ->>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
346 314 std::cout << "(srcEng=" << srcEng;
347 315 std::cout << ", srcTime=" << srcTime << "):";
348 316 std::cout << " NaN/Inf encountered";
349 317 std::cout << " (value=" << value;
350   -<<<<<<< HEAD
351 318 std::cout << ", norm=" << m_last_norm;
352 319 std::cout << ", mean=" << m_last_mean;
353 320 std::cout << ", Sigma=" << m_last_sigma;
354   -=======
355   - std::cout << ", e_norm=" << m_last_e_norm;
356   - std::cout << ", e_cut=" << m_last_e_cut;
357   - std::cout << ", power=" << m_last_power;
358   ->>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
359 321 std::cout << ")" << std::endl;
360 322 }
361 323 #endif
... ... @@ -385,6 +347,7 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
385 347 * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
386 348 * The integration is done analytically.
387 349 ***************************************************************************/
  350 +
388 351 double GModelSpectralGauss::flux(const GEnergy& emin,
389 352 const GEnergy& emax) const
390 353 {
... ... @@ -525,7 +488,6 @@ GEnergy GModelSpectralGauss::mc(const GEnergy&amp; emin,
525 488 return GEnergy(energy, "MeV");
526 489 }
527 490  
528   -
529 491 /***********************************************************************//**
530 492 * @brief Read model from XML element
531 493 *
... ... @@ -595,7 +557,6 @@ void GModelSpectralGauss::read(const GXmlElement&amp; xml)
595 557 return;
596 558 }
597 559  
598   -
599 560 /***********************************************************************//**
600 561 * @brief Write model into XML element
601 562 *
... ... @@ -834,16 +795,11 @@ double GModelSpectralGauss::eflux_kernel::eval(const double&amp; energy)
834 795 *
835 796 * Updates the precomputation cache for eval() and eval_gradients() methods.
836 797 *
837   -<<<<<<< HEAD
838   -=======
839   - * TODO integrate this to correct form
840   ->>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
841 798 ***************************************************************************/
842 799 void GModelSpectralGauss::update_eval_cache(const GEnergy& energy) const
843 800 {
844 801 // Get parameter values (takes 3 multiplications which are difficult
845 802 // to avoid)
846   -<<<<<<< HEAD
847 803 double norm = m_norm.value();
848 804 double mean = m_mean.value();
849 805 double sigma = m_sigma.value();
... ... @@ -864,31 +820,6 @@ void GModelSpectralGauss::update_eval_cache(const GEnergy&amp; energy) const
864 820 m_last_norm = eng / m_last_norm;
865 821 m_last_mean = eng / m_last_mean;
866 822 m_last_sigma = eng/ m_last_sigma;
867   -=======
868   - double index = m_index.value();
869   - double ecut = m_ecut.value();
870   - double pivot = m_pivot.value();
871   -
872   - // If the energy or one of the parameters index, cut-off or pivot
873   - // energy has changed then recompute the cache
874   - if ((m_last_energy != energy) ||
875   - (m_last_index != index) ||
876   - (m_last_ecut != ecut) ||
877   - (m_last_pivot != pivot)) {
878   -
879   - // Store actual energy and parameter values
880   - m_last_energy = energy;
881   - m_last_index = index;
882   - m_last_ecut = ecut;
883   - m_last_pivot = pivot;
884   -
885   - // Compute and store value
886   - double eng = energy.MeV();
887   - m_last_e_norm = eng / m_last_pivot;
888   - m_last_e_cut = eng / m_last_ecut;
889   - m_last_power = std::pow(m_last_e_norm, m_last_index) *
890   - std::exp(-m_last_e_cut);
891   ->>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
892 823  
893 824 } // endif: recomputation was required
894 825  
... ...