Justine Devin / gammalib

Commit 64618701770f9e297a56452d9758151c0e764eab

Authored by Jürgen Knödlseder
1 parent 647c4d11

Implement caching in GModelSpatialDiffuseMap::set_mc_cone().

Inline Side-by-side
Showing 4 changed files with 144 additions and 128 deletions
ChangeLog
  View file @6461870
1   -2016-02-16 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-02-17 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 =========================
5 5  
  6 + Optimise GModelSpatialDiffuseMap::set_mc_cone() method (#1673)
6 7 Optimise GModelSpatialDiffuseCube::set_mc_cone() method (#1673)
7 8 Use default GTime in GObservation::deadc() method (#1670)
8 9 The GApplicationPar::filename() methods now take and return GFilename (#1671)
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -16 February 2016
  5 +17 February 2016
6 6  
7 7  
8 8 1. Introduction
... ... @@ -111,10 +111,10 @@ None
111 111  
112 112 8. Model module
113 113 ---------------
114   -The GModelSpatialDiffuseCube::set_mc_cone() method now recomputes only
115   -the Monte Carlo acceptance cone in case that the cone centre or radius
116   -have changed since the last computation. This speeds up the computations
117   -(#1673).
  114 +The set_mc_cone() methods of the GModelSpatialDiffuseCube class and the
  115 +GModelSpatialDiffuseMap class now recompute only the Monte Carlo acceptance
  116 +cone in case that the cone centre or radius have changed since the last
  117 +computation. This speeds up the computations (#1673).
118 118  
119 119  
120 120 9. Numerics module
... ...
include/GModelSpatialDiffuseMap.hpp
  View file @6461870
... ... @@ -110,9 +110,11 @@ protected:
110 110 bool m_has_normalize; //!< XML has normalize attribute
111 111  
112 112 // MC simulation cache
113   - mutable double m_mc_norm; //!< Map normalization
114   - mutable std::vector<double> m_mc_cache; //!< Monte Carlo cache
115   - mutable std::vector<double> m_mc_max; //!< Monte Carlo maximum
  113 + mutable GSkyDir m_mc_centre; //!< Centre of MC cone
  114 + mutable double m_mc_radius; //!< Radius of MC cone
  115 + mutable double m_mc_norm; //!< Map normalization
  116 + mutable std::vector<double> m_mc_cache; //!< Monte Carlo cache
  117 + mutable std::vector<double> m_mc_max; //!< Monte Carlo maximum
116 118 };
117 119  
118 120  
... ...
src/model/GModelSpatialDiffuseMap.cpp
  View file @6461870
... ... @@ -654,158 +654,167 @@ void GModelSpatialDiffuseMap::write(GXmlElement&amp; xml) const
654 654 void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir& centre,
655 655 const double& radius) const
656 656 {
657   - // Initialise cache
658   - m_mc_cache.clear();
659   - m_mc_max.clear();
660   - m_mc_norm = 0.0;
  657 + // Continue only if the simulation cone has changed
  658 + if ((centre != m_mc_centre) || (radius != m_mc_radius)) {
661 659  
662   - // Determine number of skymap pixels
663   - int npix = m_map.npix();
  660 + // Save simulation cone definition
  661 + m_mc_centre = centre;
  662 + m_mc_radius = radius;
664 663  
665   - // Continue only if there are pixels
666   - if (npix > 0) {
  664 + // Initialise cache
  665 + m_mc_cache.clear();
  666 + m_mc_max.clear();
  667 + m_mc_norm = 0.0;
667 668  
668   - // Reserve space for all pixels in cache
669   - m_mc_cache.reserve(npix+1);
670   - m_mc_max.reserve(npix);
  669 + // Determine number of skymap pixels
  670 + int npix = m_map.npix();
671 671  
672   - // Set first cache value to 0
673   - m_mc_cache.push_back(0.0);
  672 + // Continue only if there are pixels
  673 + if (npix > 0) {
674 674  
675   - // Initialise cache with cumulative pixel fluxes and compute
676   - // total flux in skymap for normalization. Negative pixels are
677   - // excluded from the cumulative map. Invalid pixels are also
678   - // filtered.
679   - double sum = 0.0;
680   - double sum_map = 0.0;
681   - for (int i = 0; i < npix; ++i) {
  675 + // Reserve space for all pixels in cache
  676 + m_mc_cache.reserve(npix+1);
  677 + m_mc_max.reserve(npix);
682 678  
683   - // Derive effective pixel radius from half opening angle
684   - // that corresponds to the pixel's solid angle. For security,
685   - // the radius is enhanced by 50%.
686   - double pixel_radius =
687   - std::acos(1.0 - m_map.solidangle(i)/gammalib::twopi) *
688   - gammalib::rad2deg * 1.5;
689   -
690   - // Add up flux with simulation cone radius + effective pixel
691   - // radius. The effective pixel radius is added to make sure
692   - // that all pixels that overlap with the simulation cone are
693   - // taken into account. There is no problem of having even
694   - // pixels outside the simulation cone taken into account as
695   - // long as the mc() method has an explicit test of whether a
696   - // simulated event is contained in the simulation cone.
697   - double flux = m_map.flux(i);
698   - if (flux > 0.0) {
699   - double distance = centre.dist_deg(m_map.pix2dir(i));
700   - if (distance <= radius+pixel_radius) {
701   - sum += flux;
  679 + // Set first cache value to 0
  680 + m_mc_cache.push_back(0.0);
  681 +
  682 + // Initialise cache with cumulative pixel fluxes and compute
  683 + // total flux in skymap for normalization. Negative pixels are
  684 + // excluded from the cumulative map. Invalid pixels are also
  685 + // filtered.
  686 + double sum = 0.0;
  687 + double sum_map = 0.0;
  688 + for (int i = 0; i < npix; ++i) {
  689 +
  690 + // Derive effective pixel radius from half opening angle
  691 + // that corresponds to the pixel's solid angle. For security,
  692 + // the radius is enhanced by 50%.
  693 + double pixel_radius =
  694 + std::acos(1.0 - m_map.solidangle(i)/gammalib::twopi) *
  695 + gammalib::rad2deg * 1.5;
  696 +
  697 + // Add up flux with simulation cone radius + effective pixel
  698 + // radius. The effective pixel radius is added to make sure
  699 + // that all pixels that overlap with the simulation cone are
  700 + // taken into account. There is no problem of having even
  701 + // pixels outside the simulation cone taken into account as
  702 + // long as the mc() method has an explicit test of whether a
  703 + // simulated event is contained in the simulation cone.
  704 + double flux = m_map.flux(i);
  705 + if (flux > 0.0) {
  706 + double distance = centre.dist_deg(m_map.pix2dir(i));
  707 + if (distance <= radius+pixel_radius) {
  708 + sum += flux;
  709 + }
  710 + sum_map += flux; // sum up total flux in map
702 711 }
703   - sum_map += flux; // sum up total flux in map
704   - }
705 712  
706   - // Push back flux
707   - m_mc_cache.push_back(sum); // units: ph/cm2/s/MeV
  713 + // Push back flux
  714 + m_mc_cache.push_back(sum); // units: ph/cm2/s/MeV
708 715  
709   - } // endfor: looped over pixels
  716 + } // endfor: looped over pixels
710 717  
711   - // Normalize fluxes in the cache so that the values in the cache
712   - // run from 0 to 1.
713   - if (sum > 0.0) {
714   - for (int i = 0; i < npix; ++i) {
715   - m_mc_cache[i] /= sum;
  718 + // Normalize fluxes in the cache so that the values in the cache
  719 + // run from 0 to 1.
  720 + if (sum > 0.0) {
  721 + for (int i = 0; i < npix; ++i) {
  722 + m_mc_cache[i] /= sum;
  723 + }
716 724 }
717   - }
718 725  
719   - // Make sure that last pixel in the cache is >1
720   - m_mc_cache[npix] = 1.0001;
721   -
722   - // Set the normalization factor for the MC simulations. In case
723   - // that the map is normalized, this is the fraction of the flux
724   - // that is comprised within the simulation cone. For non-normalized
725   - // maps, this is simply the flux comprised within the simulation
726   - // cone.
727   - if (sum_map > 0.0) {
728   - if (m_normalize) {
729   - m_mc_norm = sum / sum_map;
730   - }
731   - else {
732   - m_mc_norm = sum;
  726 + // Make sure that last pixel in the cache is >1
  727 + m_mc_cache[npix] = 1.0001;
  728 +
  729 + // Set the normalization factor for the MC simulations. In case
  730 + // that the map is normalized, this is the fraction of the flux
  731 + // that is comprised within the simulation cone. For non-normalized
  732 + // maps, this is simply the flux comprised within the simulation
  733 + // cone.
  734 + if (sum_map > 0.0) {
  735 + if (m_normalize) {
  736 + m_mc_norm = sum / sum_map;
  737 + }
  738 + else {
  739 + m_mc_norm = sum;
  740 + }
733 741 }
734   - }
735 742  
736   - // Do we have a HealPix map?
737   - if (m_map.projection()->code() == "HPX") {
  743 + // Do we have a HealPix map?
  744 + if (m_map.projection()->code() == "HPX") {
738 745  
739   - // Get pointer on HealPix projection
740   - const GHealpix* healpix =
741   - static_cast<const GHealpix*>(m_map.projection());
  746 + // Get pointer on HealPix projection
  747 + const GHealpix* healpix =
  748 + static_cast<const GHealpix*>(m_map.projection());
742 749  
743   - // Compute maximum value that may occur from bilinear
744   - // interpolation within this pixel and push this value on the
745   - // stack. We do this by checking values of all neighbours.
746   - for (int i = 0; i < npix; ++i) {
  750 + // Compute maximum value that may occur from bilinear
  751 + // interpolation within this pixel and push this value on the
  752 + // stack. We do this by checking values of all neighbours.
  753 + for (int i = 0; i < npix; ++i) {
747 754  
748   - // Get neighbours
749   - std::vector<int> neighbours = healpix->neighbours(i);
  755 + // Get neighbours
  756 + std::vector<int> neighbours = healpix->neighbours(i);
750 757  
751   - // Loop over neighbours
752   - double max = m_map(i);
753   - for (int j = 0; j < neighbours.size(); ++j) {
754   - if (neighbours[j] != -1) {
755   - double value = m_map(neighbours[j]);
756   - if (value > max) {
757   - max = value;
  758 + // Loop over neighbours
  759 + double max = m_map(i);
  760 + for (int j = 0; j < neighbours.size(); ++j) {
  761 + if (neighbours[j] != -1) {
  762 + double value = m_map(neighbours[j]);
  763 + if (value > max) {
  764 + max = value;
  765 + }
758 766 }
759 767 }
760   - }
761 768  
762   - // Store maximum
763   - m_mc_max.push_back(max);
  769 + // Store maximum
  770 + m_mc_max.push_back(max);
764 771  
765   - } // endfor: looped over pixels
  772 + } // endfor: looped over pixels
766 773  
767   - } // endif: Healpix projection
  774 + } // endif: Healpix projection
768 775  
769   - // ... no, then we have a WCS map
770   - else {
  776 + // ... no, then we have a WCS map
  777 + else {
771 778  
772   - // Compute maximum value that may occur from bilinear
773   - // interpolation within this pixel and push this value on the
774   - // stack. We do this by checking the map values at the corners
775   - // and the centre of each edge.
776   - for (int i = 0; i < npix; ++i) {
777   - GSkyPixel pixel = m_map.inx2pix(i);
778   - double max = m_map(pixel);
779   - for (int ix = -1; ix < 2; ++ix) {
780   - for (int iy = -1; iy < 2; ++iy) {
781   - if (ix != 0 && iy != 0) {
782   - GSkyPixel edge(pixel.x()+ix*0.5, pixel.y()+iy*0.5);
783   - if (m_map.contains(edge)) {
784   - GSkyDir dir = m_map.pix2dir(edge);
785   - double value = m_map(dir);
786   - if (value > max) {
787   - max = value;
  779 + // Compute maximum value that may occur from bilinear
  780 + // interpolation within this pixel and push this value on the
  781 + // stack. We do this by checking the map values at the corners
  782 + // and the centre of each edge.
  783 + for (int i = 0; i < npix; ++i) {
  784 + GSkyPixel pixel = m_map.inx2pix(i);
  785 + double max = m_map(pixel);
  786 + for (int ix = -1; ix < 2; ++ix) {
  787 + for (int iy = -1; iy < 2; ++iy) {
  788 + if (ix != 0 && iy != 0) {
  789 + GSkyPixel edge(pixel.x()+ix*0.5, pixel.y()+iy*0.5);
  790 + if (m_map.contains(edge)) {
  791 + GSkyDir dir = m_map.pix2dir(edge);
  792 + double value = m_map(dir);
  793 + if (value > max) {
  794 + max = value;
  795 + }
788 796 }
789 797 }
790 798 }
791 799 }
  800 + m_mc_max.push_back(max);
792 801 }
793   - m_mc_max.push_back(max);
794   - }
795 802  
796   - } // endelse: computed maximum pixel value
  803 + } // endelse: computed maximum pixel value
797 804  
798   - // Dump cache values for debugging
799   - #if defined(G_DEBUG_CACHE)
800   - std::cout << "GModelSpatialDiffuseMap::set_mc_cone: cache";
801   - std::cout << std::endl;
802   - for (int i = 0; i < m_mc_cache.size(); ++i) {
803   - std::cout << "i=" << i;
804   - std::cout << " c=" << m_mc_cache[i] << std::endl;
805   - }
806   - #endif
  805 + // Dump cache values for debugging
  806 + #if defined(G_DEBUG_CACHE)
  807 + std::cout << "GModelSpatialDiffuseMap::set_mc_cone: cache";
  808 + std::cout << std::endl;
  809 + for (int i = 0; i < m_mc_cache.size(); ++i) {
  810 + std::cout << "i=" << i;
  811 + std::cout << " c=" << m_mc_cache[i] << std::endl;
  812 + }
  813 + #endif
  814 +
  815 + } // endif: there were cube pixels
807 816  
808   - } // endif: there were cube pixels
  817 + } // endif: simulation cone has changed
809 818  
810 819 // Return
811 820 return;
... ... @@ -918,6 +927,8 @@ void GModelSpatialDiffuseMap::init_members(void)
918 927 m_radius = 0.0;
919 928  
920 929 // Initialise MC cache
  930 + m_mc_centre.clear();
  931 + m_mc_radius = -1.0; //!< Signal initialisation
921 932 m_mc_norm = 0.0;
922 933 m_mc_cache.clear();
923 934 m_mc_max.clear();
... ... @@ -944,6 +955,8 @@ void GModelSpatialDiffuseMap::copy_members(const GModelSpatialDiffuseMap&amp; model)
944 955 m_radius = model.m_radius;
945 956  
946 957 // Copy MC cache
  958 + m_mc_centre = model.m_mc_centre;
  959 + m_mc_radius = model.m_mc_radius;
947 960 m_mc_norm = model.m_mc_norm;
948 961 m_mc_cache = model.m_mc_cache;
949 962 m_mc_max = model.m_mc_max;
... ...