Dmitry Malyshev / gammalib

Commit 94cf80555a1a61759641bf71bfa110cd98c0b47c

Authored by ellisowen
2 parents 8b56ffdf efa5ea56

Corrects the mc loop in GModelSpectralGauss.cpp, in GModelSpectralGauss::mc, as suggested

Inline Side-by-side
Showing 1 changed file with 63 additions and 0 deletions
src/model/GModelSpectralGauss.cpp
  View file @94cf805
... ... @@ -276,7 +276,11 @@ double GModelSpectralGauss::eval(const GEnergy& srcEng,
276 276 * This method simply calls the eval() method as no analytical gradients will
277 277 * be computed. See the eval() method for details.
278 278 *
  279 +<<<<<<< HEAD
279 280 * TODO: update docstring.
  281 +=======
  282 + * TODO: Implement this
  283 +>>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
280 284 ***************************************************************************/
281 285 double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng,
282 286 const GTime& srcTime)
... ... @@ -289,6 +293,7 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
289 293 // Update the evaluation cache
290 294 update_eval_cache(srcEng);
291 295  
  296 +<<<<<<< HEAD
292 297 // Compute function terms
293 298 double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
294 299 double term2 = (1 / sigma) * gammalib::inv_sqrt2pi;
... ... @@ -308,6 +313,24 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
308 313 double g_sigma = - term5 * std::exp(- term3) * (1 - (2 * term3));
309 314  
310 315 // 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
311 334 m_norm.factor_gradient(g_norm);
312 335 m_mean.factor_gradient(g_mean);
313 336 m_sigma.factor_gradient(g_sigma);
... ... @@ -315,14 +338,24 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
315 338 // Compile option: Check for NaN/Inf
316 339 #if defined(G_NAN_CHECK)
317 340 if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  341 +<<<<<<< HEAD
318 342 std::cout << "*** ERROR: GModelSpectralGauss::eval_gradients";
  343 +=======
  344 + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients";
  345 +>>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
319 346 std::cout << "(srcEng=" << srcEng;
320 347 std::cout << ", srcTime=" << srcTime << "):";
321 348 std::cout << " NaN/Inf encountered";
322 349 std::cout << " (value=" << value;
  350 +<<<<<<< HEAD
323 351 std::cout << ", norm=" << m_last_norm;
324 352 std::cout << ", mean=" << m_last_mean;
325 353 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
326 359 std::cout << ")" << std::endl;
327 360 }
328 361 #endif
... ... @@ -801,11 +834,16 @@ double GModelSpectralGauss::eflux_kernel::eval(const double&amp; energy)
801 834 *
802 835 * Updates the precomputation cache for eval() and eval_gradients() methods.
803 836 *
  837 +<<<<<<< HEAD
  838 +=======
  839 + * TODO integrate this to correct form
  840 +>>>>>>> efa5ea566f0bf64aca32f7e22a66761cbccf2f80
804 841 ***************************************************************************/
805 842 void GModelSpectralGauss::update_eval_cache(const GEnergy& energy) const
806 843 {
807 844 // Get parameter values (takes 3 multiplications which are difficult
808 845 // to avoid)
  846 +<<<<<<< HEAD
809 847 double norm = m_norm.value();
810 848 double mean = m_mean.value();
811 849 double sigma = m_sigma.value();
... ... @@ -826,6 +864,31 @@ void GModelSpectralGauss::update_eval_cache(const GEnergy&amp; energy) const
826 864 m_last_norm = eng / m_last_norm;
827 865 m_last_mean = eng / m_last_mean;
828 866 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
829 892  
830 893 } // endif: recomputation was required
831 894  
... ...