From 647c4d1108d64cfa1dfa37723cbc0f61a1f2184a Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 17 Feb 2016 00:14:40 +0100 Subject: [PATCH] GModelSpatialDiffuseCube::set_mc_cone() recomputes MC cache only in case of cone changes. --- ChangeLog | 1 + NEWS | 5 ++++- include/GModelSpatialDiffuseCube.hpp | 2 ++ src/model/GModelSpatialDiffuseCube.cpp | 310 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------- 4 files changed, 170 insertions(+), 148 deletions(-) diff --git a/ChangeLog b/ChangeLog index d1ce26d..30de987 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.1.0 released ========================= + Optimise GModelSpatialDiffuseCube::set_mc_cone() method (#1673) Use default GTime in GObservation::deadc() method (#1670) The GApplicationPar::filename() methods now take and return GFilename (#1671) Add GFilename methods and remove equivalents from GTools (#1666) diff --git a/NEWS b/NEWS index c672ea9..42e26f5 100644 --- a/NEWS +++ b/NEWS @@ -111,7 +111,10 @@ None 8. Model module --------------- -None +The GModelSpatialDiffuseCube::set_mc_cone() method now recomputes only +the Monte Carlo acceptance cone in case that the cone centre or radius +have changed since the last computation. This speeds up the computations +(#1673). 9. Numerics module diff --git a/include/GModelSpatialDiffuseCube.hpp b/include/GModelSpatialDiffuseCube.hpp index a75c4e3..295d201 100644 --- a/include/GModelSpatialDiffuseCube.hpp +++ b/include/GModelSpatialDiffuseCube.hpp @@ -130,6 +130,8 @@ protected: GEbounds m_ebounds; //!< Energy bounds of the maps // Monte Carlo cache + mutable GSkyDir m_mc_centre; //!< Centre of MC cone + mutable double m_mc_radius; //!< Radius of MC cone mutable std::vector m_mc_cache; //!< Monte Carlo cache mutable std::vector m_mc_max; //!< Maximum values for MC mutable GModelSpectralNodes m_mc_spectrum; //!< Map cube spectrum diff --git a/src/model/GModelSpatialDiffuseCube.cpp b/src/model/GModelSpatialDiffuseCube.cpp index f500ab6..6c6b74e 100644 --- a/src/model/GModelSpatialDiffuseCube.cpp +++ b/src/model/GModelSpatialDiffuseCube.cpp @@ -782,186 +782,198 @@ GEnergies GModelSpatialDiffuseCube::energies(void) void GModelSpatialDiffuseCube::set_mc_cone(const GSkyDir& centre, const double& radius) const { - // Initialise cache - m_mc_cache.clear(); - m_mc_max.clear(); - m_mc_spectrum.clear(); + // Continue only if the simulation cone has changed + if ((centre != m_mc_centre) || (radius != m_mc_radius)) { - // Fetch cube - fetch_cube(); + // Save simulation cone definition + m_mc_centre = centre; + m_mc_radius = radius; - // Determine number of cube pixels and maps - int npix = pixels(); - int nmaps = maps(); - - // Continue only if there are pixels and maps - if (npix > 0 && nmaps > 0) { - - // Reserve space for all pixels in cache - m_mc_cache.reserve((npix+1)*nmaps); - m_mc_max.reserve((npix+1)*nmaps); - m_mc_spectrum.reserve(nmaps); - - // Loop over all maps - for (int i = 0; i < nmaps; ++i) { - - // Compute pixel offset - int offset = i * (npix+1); - - // Set first cache value to 0 - m_mc_cache.push_back(0.0); - - // Initialise cache with cumulative pixel fluxes and compute - // total flux in skymap for normalization. Negative pixels are - // excluded from the cumulative map. - double sum = 0.0; - double sum_map = 0.0; - for (int k = 0; k < npix; ++k) { - - // Derive effective pixel radius from half opening angle - // that corresponds to the pixel's solid angle. For security, - // the radius is enhanced by 50%. - /* - double pixel_radius = - std::acos(1.0 - m_cube.solidangle(k)/gammalib::twopi) * - gammalib::rad2deg * 1.5; - */ - - // Add up flux with simulation cone radius + effective pixel - // radius. The effective pixel radius is added to make sure - // that all pixels that overlap with the simulation cone are - // taken into account. There is no problem of having even - // pixels outside the simulation cone taken into account as - // long as the mc() method has an explicit test of whether a - // simulated event is contained in the simulation cone. - double flux = m_cube.flux(k,i); - if (flux > 0.0) { - //double distance = centre.dist_deg(m_cube.pix2dir(k)); - //if (distance <= radius+pixel_radius) { - sum += flux; - //} - sum_map += flux; // sum up total flux in map - } + // Initialise cache + m_mc_cache.clear(); + m_mc_max.clear(); + m_mc_spectrum.clear(); - // Push back flux - m_mc_cache.push_back(sum); // units: ph/cm2/s/MeV - } + // Fetch cube + fetch_cube(); - // Normalize cumulative pixel fluxes so that the values in the - // cache run from 0 to 1 - if (sum > 0.0) { - for (int k = 0; k < npix; ++k) { - m_mc_cache[k+offset] /= sum; - } - } + // Determine number of cube pixels and maps + int npix = pixels(); + int nmaps = maps(); + + // Continue only if there are pixels and maps + if (npix > 0 && nmaps > 0) { + + // Reserve space for all pixels in cache + m_mc_cache.reserve((npix+1)*nmaps); + m_mc_max.reserve((npix+1)*nmaps); + m_mc_spectrum.reserve(nmaps); + + // Loop over all maps + for (int i = 0; i < nmaps; ++i) { - // Make sure that last pixel in the cache is >1 - m_mc_cache[npix+offset] = 1.0001; + // Compute pixel offset + int offset = i * (npix+1); - // Store flux as spectral node - if (m_logE.size() == nmaps) { + // Set first cache value to 0 + m_mc_cache.push_back(0.0); - // Compute mean energy - GEnergy energy; - energy.log10MeV(m_logE[i]); + // Initialise cache with cumulative pixel fluxes and compute + // total flux in skymap for normalization. Negative pixels + // are excluded from the cumulative map. + double sum = 0.0; + double sum_map = 0.0; + for (int k = 0; k < npix; ++k) { + + // Derive effective pixel radius from half opening angle + // that corresponds to the pixel's solid angle. For + // security, the radius is enhanced by 50%. + /* + double pixel_radius = + std::acos(1.0 - m_cube.solidangle(k)/gammalib::twopi) * + gammalib::rad2deg * 1.5; + */ + + // Add up flux with simulation cone radius + effective + // pixel radius. The effective pixel radius is added to + // make sure that all pixels that overlap with the + // simulation cone are taken into account. There is no + // problem of having even pixels outside the simulation + // cone taken into account as long as the mc() method has + // an explicit test of whether a simulated event is + // contained in the simulation cone. + double flux = m_cube.flux(k,i); + if (flux > 0.0) { + //double distance = centre.dist_deg(m_cube.pix2dir(k)); + //if (distance <= radius+pixel_radius) { + sum += flux; + //} + sum_map += flux; // sum up total flux in map + } - // Only append node if flux > 0 + // Push back flux + m_mc_cache.push_back(sum); // units: ph/cm2/s/MeV + } + + // Normalize cumulative pixel fluxes so that the values in + // the cache run from 0 to 1 if (sum > 0.0) { - m_mc_spectrum.append(energy, sum); + for (int k = 0; k < npix; ++k) { + m_mc_cache[k+offset] /= sum; + } } - } + // Make sure that last pixel in the cache is >1 + m_mc_cache[npix+offset] = 1.0001; - // Do we have a HealPix map? - if (m_cube.projection()->code() == "HPX") { + // Store flux as spectral node + if (m_logE.size() == nmaps) { - // Get pointer on HealPix projection - const GHealpix* healpix = - static_cast(m_cube.projection()); + // Compute mean energy + GEnergy energy; + energy.log10MeV(m_logE[i]); - // Compute maximum value that may occur from bilinear - // interpolation within this pixel and push this value on the - // stack. We do this by checking values of all neighbours. - for (int k = 0; k < npix; ++k) { + // Only append node if flux > 0 + if (sum > 0.0) { + m_mc_spectrum.append(energy, sum); + } - // Get neighbours - std::vector neighbours = healpix->neighbours(k); + } - // Loop over neighbours - double max = m_cube(k,i); - for (int j = 0; j < neighbours.size(); ++j) { - if (neighbours[j] != -1) { - double value = m_cube(neighbours[j],i); - if (value > max) { - max = value; + // Do we have a HealPix map? + if (m_cube.projection()->code() == "HPX") { + + // Get pointer on HealPix projection + const GHealpix* healpix = + static_cast(m_cube.projection()); + + // Compute maximum value that may occur from bilinear + // interpolation within this pixel and push this value + // on the stack. We do this by checking values of all + // neighbours. + for (int k = 0; k < npix; ++k) { + + // Get neighbours + std::vector neighbours = healpix->neighbours(k); + + // Loop over neighbours + double max = m_cube(k,i); + for (int j = 0; j < neighbours.size(); ++j) { + if (neighbours[j] != -1) { + double value = m_cube(neighbours[j],i); + if (value > max) { + max = value; + } } } - } - // Store maximum - m_mc_max[k+offset] = max; + // Store maximum + m_mc_max[k+offset] = max; - } // endfor: looped over pixels + } // endfor: looped over pixels - } // endif: Healpix projection + } // endif: Healpix projection - // ... no, then we have a WCS map - else { + // ... no, then we have a WCS map + else { - // Compute maximum value that may occur from bilinear - // interpolation within this pixel and push this value on the - // stack. We do this by checking the map values at the corners - // and the centre of each edge. - for (int k = 0; k < npix; ++k) { - GSkyPixel pixel = m_cube.inx2pix(k); - double max = m_cube(pixel,i); - for (int ix = -1; ix < 2; ++ix) { - for (int iy = -1; iy < 2; ++iy) { - if (ix != 0 && iy != 0) { - GSkyPixel edge(pixel.x()+ix*0.5, pixel.y()+iy*0.5); - if (m_cube.contains(edge)) { - GSkyDir dir = m_cube.pix2dir(edge); - double value = m_cube(dir,i); - if (value > max) { - max = value; + // Compute maximum value that may occur from bilinear + // interpolation within this pixel and push this value + // on the stack. We do this by checking the map values + // at the corners and the centre of each edge. + for (int k = 0; k < npix; ++k) { + GSkyPixel pixel = m_cube.inx2pix(k); + double max = m_cube(pixel,i); + for (int ix = -1; ix < 2; ++ix) { + for (int iy = -1; iy < 2; ++iy) { + if (ix != 0 && iy != 0) { + GSkyPixel edge(pixel.x()+ix*0.5, pixel.y()+iy*0.5); + if (m_cube.contains(edge)) { + GSkyDir dir = m_cube.pix2dir(edge); + double value = m_cube(dir,i); + if (value > max) { + max = value; + } } } } } + m_mc_max[k+offset] = max; } - m_mc_max[k+offset] = max; - } - } // endelse: computed maximum pixel value + } // endelse: computed maximum pixel value - // Make sure that last pixel in the maximum cache has a well defined value - m_mc_max[npix+offset] = 0.0; + // Make sure that last pixel in the maximum cache has a + // well defined value + m_mc_max[npix+offset] = 0.0; - } // endfor: looped over all maps + } // endfor: looped over all maps - // Dump cache values for debugging - #if defined(G_DEBUG_CACHE) - std::cout << "GModelSpatialDiffuseCube::set_mc_cone: cache"; - std::cout << std::endl; - for (int i = 0; i < m_mc_cache.size(); ++i) { - std::cout << "i=" << i; - std::cout << " cache=" << m_mc_cache[i]; - std::cout << " max=" << m_mc_max[i] << std::endl; - } - #endif - - // Dump spectrum for debugging - #if defined(G_DEBUG_CACHE) - std::cout << "GModelSpatialDiffuseCube::set_mc_cone: spectrum"; - std::cout << std::endl; - for (int i = 0; i < m_mc_spectrum.nodes(); ++i) { - std::cout << i; - std::cout << " " << m_mc_spectrum.energy(i); - std::cout << " " << m_mc_spectrum.intensity(i) << std::endl; - } - #endif + // Dump cache values for debugging + #if defined(G_DEBUG_CACHE) + std::cout << "GModelSpatialDiffuseCube::set_mc_cone: cache"; + std::cout << std::endl; + for (int i = 0; i < m_mc_cache.size(); ++i) { + std::cout << "i=" << i; + std::cout << " cache=" << m_mc_cache[i]; + std::cout << " max=" << m_mc_max[i] << std::endl; + } + #endif + + // Dump spectrum for debugging + #if defined(G_DEBUG_CACHE) + std::cout << "GModelSpatialDiffuseCube::set_mc_cone: spectrum"; + std::cout << std::endl; + for (int i = 0; i < m_mc_spectrum.nodes(); ++i) { + std::cout << i; + std::cout << " " << m_mc_spectrum.energy(i); + std::cout << " " << m_mc_spectrum.intensity(i) << std::endl; + } + #endif + + } // endif: there were cube pixels and maps - } // endif: there were cube pixels and maps + } // endif: simulation cone has changed // Return return; @@ -1155,6 +1167,8 @@ void GModelSpatialDiffuseCube::init_members(void) m_loaded = false; // Initialise MC cache + m_mc_centre.clear(); + m_mc_radius = -1.0; //!< Signal that initialisation is needed m_mc_cache.clear(); m_mc_max.clear(); m_mc_spectrum.clear(); @@ -1180,6 +1194,8 @@ void GModelSpatialDiffuseCube::copy_members(const GModelSpatialDiffuseCube& mode m_loaded = model.m_loaded; // Copy MC cache + m_mc_centre = model.m_mc_centre; + m_mc_radius = model.m_mc_radius; m_mc_cache = model.m_mc_cache; m_mc_max = model.m_mc_max; m_mc_spectrum = model.m_mc_spectrum; -- libgit2 0.21.4