diff --git a/include/GSkyMap.hpp b/include/GSkyMap.hpp index e848293..44caac9 100644 --- a/include/GSkyMap.hpp +++ b/include/GSkyMap.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GSkyMap.hpp - Sky map class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2016 by Juergen Knoedlseder * + * copyright (C) 2010-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -85,6 +85,10 @@ class GVector; class GSkyMap : public GBase { friend GSkyMap sqrt(const GSkyMap& map); + friend GSkyMap log(const GSkyMap& map); + friend GSkyMap log10(const GSkyMap& map); + friend GSkyMap abs(const GSkyMap& map); + friend GSkyMap sign(const GSkyMap& map); public: // Constructors and destructors diff --git a/pyext/GSkyMap.i b/pyext/GSkyMap.i index 2ed6ae1..329f0a2 100644 --- a/pyext/GSkyMap.i +++ b/pyext/GSkyMap.i @@ -245,6 +245,18 @@ public: GSkyMap sqrt() { return sqrt(*self); } + GSkyMap log() { + return log(*self); + } + GSkyMap log10() { + return log10(*self); + } + GSkyMap abs() { + return abs(*self); + } + GSkyMap sign() { + return sign(*self); + } // Python 2.x operator/= GSkyMap __idiv__(const GSkyMap& map) { self->operator/=(map); diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index 12d477d..3b22fe3 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -76,6 +76,9 @@ #define G_READ_HEALPIX "GSkyMap::read_healpix(GFitsTable*)" #define G_READ_WCS "GSkyMap::read_wcs(GFitsImage*)" #define G_ALLOC_WCS "GSkyMap::alloc_wcs(GFitsImage*)" +#define G_SQRT "GSkyMap sqrt(GSkyMap&)" +#define G_LOG "GSkyMap log(GSkyMap&)" +#define G_LOG10 "GSkyMap log10(GSkyMap&)" /* __ Macros _____________________________________________________________ */ @@ -3019,10 +3022,10 @@ bool GSkyMap::is_wcs(const GFitsHDU& hdu) const ==========================================================================*/ /***********************************************************************//** - * @brief Computes square root of skymap elements + * @brief Computes square root of sky map elements * - * @param[in] map Skymap. - * @return Skymap containing the square root of every element. + * @param[in] map Sky map. + * @return Sky map containing the square root of every element. ***************************************************************************/ GSkyMap sqrt(const GSkyMap& map) { @@ -3038,20 +3041,177 @@ GSkyMap sqrt(const GSkyMap& map) // Get the content from the bin double content = map(j,i); - // Check if content is not negative + // Throw an exception if content is negative if (content < 0.0) { - std::string msg = "Negative value envountered." - " Cannot take the sqrt from a negative value"; - throw GException::invalid_value(G_OP_UNARY_DIV, msg); + std::string msg = "Negative value encountered. " + "Cannot take the sqrt from a negative " + "value"; + throw GException::invalid_value(G_SQRT, msg); } // Set content of the result map result(j,i) = std::sqrt(content); - } // endfor: Loop over all bins + } // endfor: looped over all bins + + } // endfor: looped over all maps + + // Return sky map + return result; +} + + +/***********************************************************************//** + * @brief Computes the natural logarithm of sky map elements + * + * @param[in] map Sky map. + * @return Sky map containing the natural logarithm of every element. + ***************************************************************************/ +GSkyMap log(const GSkyMap& map) +{ + // Initialise result vector + GSkyMap result(map); + + // Loop over all maps + for (int i = 0; i < map.nmaps(); ++i) { + + // Loop over all bins + for (int j = 0; j < map.npix(); ++j) { + + // Get the content from the bin + double content = map(j,i); + + // Throw an exception if content is negative + if (content <= 0.0) { + std::string msg = "Negative or zero value encountered. " + "Cannot calculate the logarithm of a " + "negative or zero value"; + throw GException::invalid_value(G_LOG, msg); + } + + // Set content of the result map + result(j,i) = std::log(content); + + } // endfor: looped over all bins + + } // endfor: looped over all maps + + // Return sky map + return result; +} + + +/***********************************************************************//** + * @brief Computes the base 10 logarithm of sky map elements + * + * @param[in] map Sky map. + * @return Sky map containing the base 10 logarithm of every element. + ***************************************************************************/ +GSkyMap log10(const GSkyMap& map) +{ + // Initialise result vector + GSkyMap result(map); + + // Loop over all maps + for (int i = 0; i < map.nmaps(); ++i) { + + // Loop over all bins + for (int j = 0; j < map.npix(); ++j) { + + // Get the content from the bin + double content = map(j,i); + + // Throw an exception if content is negative + if (content <= 0.0) { + std::string msg = "Negative or zero value encountered. " + "Cannot calculate the logarithm of a " + "negative or zero value"; + throw GException::invalid_value(G_LOG10, msg); + } + + // Set content of the result map + result(j,i) = std::log10(content); + + } // endfor: looped over all bins + + } // endfor: looped over all maps + + // Return sky map + return result; +} + + +/***********************************************************************//** + * @brief Computes the absolute value of sky map elements + * + * @param[in] map Sky map. + * @return Sky map containing the absolute value of every element. + ***************************************************************************/ +GSkyMap abs(const GSkyMap& map) +{ + // Initialise result vector + GSkyMap result(map); + + // Loop over all maps + for (int i = 0; i < map.nmaps(); ++i) { + + // Loop over all bins + for (int j = 0; j < map.npix(); ++j) { + + // Get the content from the bin + double content = map(j,i); + + // Set content of the result map + result(j,i) = std::abs(content); + + } // endfor: looped over all bins + + } // endfor: looped over all maps + + // Return sky map + return result; +} + + +/***********************************************************************//** + * @brief Computes the sign value of sky map elements + * + * @param[in] map Sky map. + * @return Sky map containing the sign value of every pixel. + * + * This method returns a sky map filled with a value of 1 if the pixel + * is positive, a value of -1 if the pixel is negative or a value of 0 + * if the pixel is 0. + ***************************************************************************/ +GSkyMap sign(const GSkyMap& map) +{ + // Initialise result vector + GSkyMap result(map); + + // Loop over all maps + for (int i = 0; i < map.nmaps(); ++i) { + + // Loop over all bins + for (int j = 0; j < map.npix(); ++j) { + + // Get the content from the bin + double content = map(j,i); + + // Handle the 3 cases + if (content < 0.0) { + result(j,i) = -1.0; + } + else if (content > 0.0) { + result(j,i) = +1.0; + } + else { + result(j,i) = 0.0; + } + + } // endfor: looped over all bins - } // endfor: Loop over all maps + } // endfor: looed over all maps - // Return vector + // Return sky map return result; } diff --git a/test/test_GSky.py b/test/test_GSky.py index fd14630..bde4433 100755 --- a/test/test_GSky.py +++ b/test/test_GSky.py @@ -72,6 +72,7 @@ class Test(gammalib.GPythonTestSuite): # Append tests self.append(self.test_methods, "Test sky map methods") + self.append(self.test_friend_methods, "Test sky map friend methods") self.append(self.test_operators, "Test sky map operators") self.append(self.test_skymap_healpix, "Test HEALPix map") self.append(self.test_skymap_ait, "Test AIT projection map") @@ -123,7 +124,35 @@ class Test(gammalib.GPythonTestSuite): # Return return - + + # Test sky map friend methods + def test_friend_methods(self): + """ + Test the sky map friend methods + """ + # Setup sky maps + m = gammalib.GSkyMap('CAR', 'CEL', 83.63, 22.01, -3.7, 2.6, 10, 8, 12) + + # Fill all values of skymaps with numbers > 0 + values = range(1,m.npix()+1) + for j in range(m.nmaps()): + for i in range(m.npix()): + m[i,j] = values[i] + + # Test map values after calling friend methods on one pixel + ipix = 18 + self.test_value(m.sqrt()[ipix,0], values[ipix]**0.5, 'Check sqrt method') + self.test_value(m.log()[ipix,0], math.log(values[ipix]), 'Check log method') + self.test_value(m.log10()[ipix,0], math.log10(values[ipix]), 'Check log10 method') + + # Set the interesting value to a negative value + m[ipix,0] = -values[ipix] + self.test_value(m.abs()[ipix,0], values[ipix], 'Check abs method') + self.test_value(m.sign()[ipix,0], -1, 'Check sign method') + + # Return + return + # Test sky map operators def test_operators(self): """