Justine Devin / gammalib

Commit 90b77b79e764fd19cc5d1788d2f77eefc381c0cd

Authored by Michael Mayer
1 parent 8a0280de

implemented isconstant()-method and caching logic for eval()

Inline Side-by-side
Showing 4 changed files with 183 additions and 56 deletions
inst/cta/include/GCTAModelBackground.hpp
  View file @90b77b7
... ... @@ -67,6 +67,7 @@ public:
67 67 virtual void clear(void);
68 68 virtual GCTAModelBackground* clone(void) const;
69 69 virtual std::string type(void) const { return "CTABackground"; }
  70 + virtual bool isconstant(void) const;
70 71 virtual double eval(const GEvent& event,
71 72 const GObservation& obs) const;
72 73 virtual double eval_gradients(const GEvent& event,
... ... @@ -158,6 +159,27 @@ protected:
158 159 GModelSpectral* m_spectral; //!< Spectral model
159 160 GModelTemporal* m_temporal; //!< Temporal model
160 161 GMatrix m_rot; //!!< Rotation matrix from model system to skydir
  162 +
  163 + // Npred cache
  164 + mutable std::vector<std::string> m_npred_names; //!< Model names
  165 + mutable std::vector<GEnergy> m_npred_energies; //!< Model energy
  166 + mutable std::vector<GTime> m_npred_times; //!< Model time
  167 + mutable std::vector<double> m_npred_values; //!< Model values
  168 +
161 169 };
162 170  
163 171 #endif /* GMODELSPATIAL_HPP */
  172 +
  173 +/***********************************************************************//**
  174 + * @brief Signals if sky model is temporally constant
  175 + *
  176 + * @return True if sky model is temporally constant, false otherwise.
  177 + *
  178 + * Signals if the sky model is temporally constant. A temporally constant
  179 + * model is a model that has a temporal component of type "Constant".
  180 + ***************************************************************************/
  181 +inline
  182 +bool GCTAModelBackground::isconstant(void) const
  183 +{
  184 + return (m_temporal != NULL && m_temporal->type() == "Constant");
  185 +}
... ...
inst/cta/pyext/GCTAModelBackground.i
  View file @90b77b7
... ... @@ -50,6 +50,7 @@ public:
50 50 virtual void clear(void);
51 51 virtual GCTAModelBackground* clone(void) const;
52 52 virtual std::string type(void) const;
  53 + virtual bool isconstant(void) const;
53 54 virtual double eval(const GEvent& event,
54 55 const GObservation& obs) const;
55 56 virtual double eval_gradients(const GEvent& event,
... ...
inst/cta/src/GCTAModelBackground.cpp
  View file @90b77b7
... ... @@ -69,8 +69,10 @@ const GModelRegistry g_cta_model_background_registry(&amp;g_cta_model_bac
69 69 /* __ Macros _____________________________________________________________ */
70 70  
71 71 /* __ Coding definitions _________________________________________________ */
  72 +#define G_USE_NPRED_CACHE //!< Use Npred cache in npred_diffuse method
72 73  
73 74 /* __ Debug definitions __________________________________________________ */
  75 +//#define G_DEBUG_NPRED //!< Debug npred_diffuse method
74 76 //#define G_DUMP_MC //!< Dump MC information
75 77  
76 78  
... ... @@ -449,67 +451,123 @@ double GCTAModelBackground::npred(const GEnergy&amp; obsEng,
449 451 {
450 452 // Initialise result
451 453 double npred = 0.0;
  454 + bool has_npred = false;
  455 +
  456 + // Build unique identifier
  457 + std::string id = obs.instrument() + "::" + obs.id();
  458 + //
  459 + // Check if Npred value is already in cache
  460 + #if defined(G_USE_NPRED_CACHE)
  461 + if (!m_npred_names.empty()) {
  462 +
  463 + // Search for unique identifier, and if found, recover Npred value
  464 + // and break
  465 + for (int i = 0; i < m_npred_names.size(); ++i) {
  466 + if (m_npred_names[i] == id &&
  467 + m_npred_energies[i] == obsEng) {
  468 + npred = m_npred_values[i];
  469 + has_npred = true;
  470 + #if defined(G_DEBUG_NPRED)
  471 + std::cout << "GCTAModelBackground::npred:";
  472 + std::cout << " cache=" << i;
  473 + std::cout << " npred=" << npred << std::endl;
  474 + #endif
  475 + break;
  476 + }
  477 + }
452 478  
453   - // Evaluate only if model is valid
454   - if (valid_model()) {
455   -
456   - // Get pointer on CTA events list
457   - const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
458   - if (events == NULL) {
459   - throw GException::no_list(G_NPRED);
460   - }
461   -
462   - //todo we need to think about if we always want to integrate
463   - // starting from the pointing direction or instead take
464   - // the center of gravity of the model
465   - // For now we use the pointing position as integration start
466   -
467   - // Get CTA pointing direction
468   - GCTAPointing* pnt = dynamic_cast<GCTAPointing*>(obs.pointing());
469   - if (pnt == NULL) {
470   - throw GCTAException::no_pointing(G_NPRED);
471   - }
  479 + } // endif: there were values in the Npred cache
  480 + #endif
  481 + //std::cout<<"\nCALLING NPRED "<<id<<" ...HAD NPRED: "<<has_npred<<std::endl;
472 482  
473   - // Get ROI radius in radians
474   - double roi_radius = events->roi().radius() * gammalib::deg2rad;
  483 + // Continue only if no Npred cache value was found
  484 + if (!has_npred) {
475 485  
476   - // Get distance from ROI centre in radians
477   - double roi_distance = events->roi().centre().dist(pnt->dir());
  486 + // Evaluate only if model is valid
  487 + if (valid_model()) {
478 488  
479   - // Initialise rotation matrix to transform from model system to celestial system
480   - GMatrix ry;
481   - GMatrix rz;
482   - ry.eulery(events->roi().centre().dec_deg() - 90.0);
483   - rz.eulerz(-events->roi().centre().ra_deg());
484   - GMatrix rot = (ry * rz).transpose();
  489 + // Get pointer on CTA events list
  490 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  491 + if (events == NULL) {
  492 + throw GException::no_list(G_NPRED);
  493 + }
485 494  
486   - // Compute position angle of ROI centre with respect to model
487   - // centre (radians)
488   - double omega0 = pnt->dir().posang(events->roi().centre().dir());
  495 + //todo we need to think about if we always want to integrate
  496 + // starting from the pointing direction or instead take
  497 + // the center of gravity of the model
  498 + // For now we use the pointing position as integration start
489 499  
490   - // Setup integration function
491   - GCTAModelBackground::npred_roi_kern_theta integrand(spatial(), obsEng,obsTime,rot, roi_radius, roi_distance,omega0);
  500 + // Get CTA pointing direction
  501 + GCTAPointing* pnt = dynamic_cast<GCTAPointing*>(obs.pointing());
  502 + if (pnt == NULL) {
  503 + throw GCTAException::no_pointing(G_NPRED);
  504 + }
492 505  
493   - // Setup integrator
494   - GIntegral integral(&integrand);
  506 + // Get ROI radius in radians
  507 + double roi_radius = events->roi().radius() * gammalib::deg2rad;
495 508  
496   - // Setup integration boundaries
497   - double rmin = (roi_distance > roi_radius) ? roi_distance-roi_radius : 0.0;
498   - double rmax = roi_radius + roi_distance;
  509 + // Get distance from ROI centre in radians
  510 + double roi_distance = events->roi().centre().dist(pnt->dir());
499 511  
500   - // Spatially integrate radial component
501   - npred = integral.romb(rmin, rmax);
  512 + // Initialise rotation matrix to transform from model system to celestial system
  513 + GMatrix ry;
  514 + GMatrix rz;
  515 + ry.eulery(events->roi().centre().dec_deg() - 90.0);
  516 + rz.eulerz(-events->roi().centre().ra_deg());
  517 + GMatrix rot = (ry * rz).transpose();
502 518  
503   - // Multiply in spectral and temporal components
504   - npred *= spectral()->eval(obsEng, obsTime);
505   - npred *= temporal()->eval(obsTime);
  519 + // Compute position angle of ROI centre with respect to model
  520 + // centre (radians)
  521 + double omega0 = pnt->dir().posang(events->roi().centre().dir());
506 522  
507   - // Apply deadtime correction
508   - npred *= obs.deadc(obsTime);
  523 + // Setup integration function
  524 + GCTAModelBackground::npred_roi_kern_theta integrand(spatial(), obsEng,obsTime,rot, roi_radius, roi_distance,omega0);
509 525  
510   - } // endif: model was valid
  526 + // Setup integrator
  527 + GIntegral integral(&integrand);
  528 + integral.eps(1e-2);
  529 +
  530 + // Setup integration boundaries
  531 + double rmin = (roi_distance > roi_radius) ? roi_distance-roi_radius : 0.0;
  532 + double rmax = roi_radius + roi_distance;
  533 +
  534 + // Spatially integrate radial component
  535 + npred = integral.romb(rmin, rmax);
  536 +
  537 + // Store result in Npred cache
  538 + #if defined(G_USE_NPRED_CACHE)
  539 + m_npred_names.push_back(id);
  540 + m_npred_energies.push_back(obsEng);
  541 + m_npred_times.push_back(obsTime);
  542 + m_npred_values.push_back(npred);
  543 + #endif
  544 +
  545 + // Debug: Check for NaN
  546 + #if defined(G_NAN_CHECK)
  547 + if (gammalib::isnotanumber(npred) || gammalib::isinfinite(npred)) {
  548 + std::cout << "*** ERROR: GCTAModelBackground::npred:";
  549 + std::cout << " NaN/Inf encountered";
  550 + std::cout << " (npred=" << npred;
  551 + std::cout << ", roi_radius=" << roi_radius;
  552 + std::cout << ")" << std::endl;
  553 + }
  554 + #endif
  555 +
  556 + } // endif: model was valid
  557 +
  558 + } // endif: Npred computation required
  559 +
  560 + //if(obsEng.TeV()>0.795&& obsEng.TeV()<1.1){std::cout<<"NPRED SPATIAL = "<<npred<<std::endl;}
  561 +
  562 + // Multiply in spectral and temporal components
  563 + npred *= spectral()->eval(obsEng, obsTime);
  564 + npred *= temporal()->eval(obsTime);
  565 + //std::cout<<npred<<" ohne deadc ";
  566 + // Apply deadtime correction
  567 + npred *= obs.deadc(obsTime);
  568 + //std::cout<<npred<<std::endl;
511 569  
512   - // Return
  570 + // Return Npred
513 571 return npred;
514 572 }
515 573  
... ... @@ -737,7 +795,7 @@ void GCTAModelBackground::write(GXmlElement&amp; xml) const
737 795 spectral()->write(*spec);
738 796 }
739 797  
740   - // Write radial model
  798 + // Write spatial model
741 799 if (spatial()) {
742 800 GXmlElement* spat = src->element("spatialModel", 0);
743 801 spatial()->write(*spat);
... ... @@ -843,6 +901,12 @@ void GCTAModelBackground::init_members(void)
843 901 m_spectral = NULL;
844 902 m_temporal = NULL;
845 903  
  904 + // Initialise Npred cache
  905 + m_npred_names.clear();
  906 + m_npred_energies.clear();
  907 + m_npred_times.clear();
  908 + m_npred_values.clear();
  909 +
846 910 // Return
847 911 return;
848 912 }
... ... @@ -857,6 +921,12 @@ void GCTAModelBackground::init_members(void)
857 921 ***************************************************************************/
858 922 void GCTAModelBackground::copy_members(const GCTAModelBackground& model)
859 923 {
  924 + // Copy cache
  925 + m_npred_names = model.m_npred_names;
  926 + m_npred_energies = model.m_npred_energies;
  927 + m_npred_times = model.m_npred_times;
  928 + m_npred_values = model.m_npred_values;
  929 +
860 930 // Clone radial, spectral and temporal model components
861 931 m_spatial = (model.m_spatial != NULL) ? model.m_spatial->clone() : NULL;
862 932 m_spectral = (model.m_spectral != NULL) ? model.m_spectral->clone() : NULL;
... ... @@ -877,6 +947,11 @@ void GCTAModelBackground::copy_members(const GCTAModelBackground&amp; model)
877 947 ***************************************************************************/
878 948 void GCTAModelBackground::free_members(void)
879 949 {
  950 + m_npred_names.clear();
  951 + m_npred_energies.clear();
  952 + m_npred_times.clear();
  953 + m_npred_values.clear();
  954 +
880 955 // Free memory
881 956 if (m_spatial != NULL) delete m_spatial;
882 957 if (m_spectral != NULL) delete m_spectral;
... ... @@ -1109,6 +1184,9 @@ double GCTAModelBackground::npred_roi_kern_theta::eval(double theta)
1109 1184  
1110 1185 // Integrate over phi
1111 1186 GIntegral integral(&integrand);
  1187 + integral.eps(1e-2);
  1188 + //std::cout<<omega_min<<" "<<omega_max<<std::endl;
  1189 + //value = integral.romb(0.0,gammalib::twopi) * sin_theta;
1112 1190 value = integral.romb(omega_min,omega_max) * sin_theta;
1113 1191  
1114 1192 // Debug: Check for NaN
... ...
src/model/GModelSpatialDiffuseCube.cpp
  View file @90b77b7
... ... @@ -54,7 +54,7 @@ const GModelSpatialRegistry g_spatial_cube_registry(&amp;g_spatial_cube_seed);
54 54 /* __ Coding definitions _________________________________________________ */
55 55  
56 56 /* __ Debug definitions __________________________________________________ */
57   -
  57 +#define G_DEBUG_PREPARE
58 58  
59 59 /*==========================================================================
60 60 = =
... ... @@ -279,8 +279,21 @@ GModelSpatialDiffuseCube* GModelSpatialDiffuseCube::clone(void) const
279 279 ***************************************************************************/
280 280 double GModelSpatialDiffuseCube::eval(const GPhoton& photon) const
281 281 {
282   - // Get skymap bin
283   - int imap = m_ebounds.index(photon.energy());
  282 + // Initialise skymap bin
  283 + int imap = -1;
  284 +
  285 + // Get skymap bin. If energy outside bounds, use the closest bin
  286 + if( photon.energy() < m_ebounds.emin() ) {
  287 + imap = 0;
  288 + } //endif: energy was below emin
  289 +
  290 + else if(photon.energy() > m_ebounds.emax()) {
  291 + imap = m_ebounds.size()-1;
  292 + } // endif: energy was above emax
  293 +
  294 + else {
  295 + imap = m_ebounds.index(photon.energy());
  296 + } // endif: found skymap bin
284 297  
285 298 // Get skymap intensity
286 299 double intensity = m_cube(photon.dir(), imap);
... ... @@ -301,8 +314,21 @@ double GModelSpatialDiffuseCube::eval(const GPhoton&amp; photon) const
301 314 ***************************************************************************/
302 315 double GModelSpatialDiffuseCube::eval_gradients(const GPhoton& photon) const
303 316 {
304   - // Get skymap bin
305   - int imap = m_ebounds.index(photon.energy());
  317 + // Initialise skymap bin
  318 + int imap = -1;
  319 +
  320 + // Get skymap bin. If energy outside bounds, use the closest bin
  321 + if( photon.energy() < m_ebounds.emin() ) {
  322 + imap = 0;
  323 + } //endif: energy was below emin
  324 +
  325 + else if(photon.energy() > m_ebounds.emax()) {
  326 + imap = m_ebounds.size()-1;
  327 + } // endif: energy was abvoe emax
  328 +
  329 + else {
  330 + imap = m_ebounds.index(photon.energy());
  331 + } // endif: found skymap bin
306 332  
307 333 // Get skymap intensity
308 334 double intensity = m_cube(photon.dir(), imap);
... ... @@ -572,7 +598,7 @@ void GModelSpatialDiffuseCube::load(const std::string&amp; filename)
572 598  
573 599 // Get the column with the name "Energy"
574 600 GFitsTableDoubleCol* ptr_energy = (GFitsTableDoubleCol*)&(*energies)["Energy"];
575   - ptr_energy->print();
  601 + //ptr_energy->print();
576 602  
577 603 // Get the unit of the energies
578 604 // Default for Fermi Galactic diffuse model
... ... @@ -617,7 +643,7 @@ void GModelSpatialDiffuseCube::load(const std::string&amp; filename)
617 643 } // endif: ENERGIES extension existed
618 644  
619 645 // Prepare cube
620   - prepare_cube();
  646 + //prepare_cube();
621 647  
622 648 // Cube is loaded
623 649 m_loaded = true;
... ...