Justine Devin / gammalib

Commit 961a8778a36b764b2a2384f0b66b19b3ee35cb00

Authored by Jürgen Knödlseder
1 parent 0ef13cad

Add scaling operators to GSkymap class.

Inline Side-by-side
Showing 7 changed files with 162 additions and 25 deletions
ChangeLog
  View file @961a877
1   -2015-02-05 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2015-02-08 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.0.0 released
4 4 ========================
... ... @@ -13,6 +13,7 @@
13 13 Add GBilinear class
14 14 Add WCS projections (MOL)
15 15 Add GXmlElement::has_attribute() method
  16 + Add GSkymap scaling operators (#1296)
16 17  
17 18  
18 19 2015-01-07 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.0.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -5 February 2015
  5 +8 February 2015
6 6  
7 7  
8 8 1. Introduction
... ... @@ -81,6 +81,9 @@ None
81 81 The following WCS projections have been added:
82 82 * GWcsMOL (Mollweide's projection)
83 83  
  84 +Scaling operators operator*=(double) and operator/=(double) have been added
  85 +to GSkymap to allow scaling of all pixel values by a given value.
  86 +
84 87  
85 88 13. Support module
86 89 ------------------
... ...
include/GSkymap.hpp
  View file @961a877
... ... @@ -107,7 +107,9 @@ public:
107 107 GSkymap& operator+=(const GSkymap& map);
108 108 GSkymap& operator-=(const GSkymap& map);
109 109 GSkymap& operator*=(const GSkymap& map);
  110 + GSkymap& operator*=(const double& factor);
110 111 GSkymap& operator/=(const GSkymap& map);
  112 + GSkymap& operator/=(const double& factor);
111 113 double& operator()(const int& index, const int& map = 0);
112 114 const double& operator()(const int& index, const int& map = 0) const;
113 115 double& operator()(const GSkyPixel& pixel, const int& map = 0);
... ...
pyext/GSkymap.i
  View file @961a877
... ... @@ -131,6 +131,7 @@ public:
131 131 GSkymap& operator+=(const GSkymap& map);
132 132 GSkymap& operator-=(const GSkymap& map);
133 133 GSkymap& operator*=(const GSkymap& map);
  134 + GSkymap& operator*=(const double& factor);
134 135  
135 136 // Methods
136 137 void clear(void);
... ... @@ -211,4 +212,14 @@ public:
211 212 self->operator/=(map);
212 213 return (*self);
213 214 }
  215 + // Python 2.x operator/=
  216 + GSkymap __idiv__(const double& factor) {
  217 + self->operator/=(factor);
  218 + return (*self);
  219 + }
  220 + // Python 3.x operator/=
  221 + GSkymap __itruediv__(const double& factor) {
  222 + self->operator/=(factor);
  223 + return (*self);
  224 + }
214 225 };
... ...
src/sky/GSkymap.cpp
  View file @961a877
... ... @@ -357,7 +357,7 @@ GSkymap&amp; GSkymap::operator+=(const GSkymap&amp; map)
357 357  
358 358 } // endfor: looped over all pixels
359 359  
360   - // Return present sky map
  360 + // Return this object
361 361 return *this;
362 362 }
363 363  
... ... @@ -408,7 +408,7 @@ GSkymap&amp; GSkymap::operator-=(const GSkymap&amp; map)
408 408  
409 409 } // endfor: looped over all pixels
410 410  
411   - // Return present sky map
  411 + // Return this object
412 412 return *this;
413 413 }
414 414  
... ... @@ -458,7 +458,31 @@ GSkymap&amp; GSkymap::operator*=(const GSkymap&amp; map)
458 458  
459 459 } // endfor: looped over all pixels
460 460  
461   - // Return present sky map
  461 + // Return this object
  462 + return *this;
  463 +}
  464 +
  465 +
  466 +/***********************************************************************//**
  467 + * @brief Sky map scaling operator
  468 + *
  469 + * @param[in] factor Scale factor.
  470 + * @return Sky map.
  471 + *
  472 + * Multiplies all pixels of the sky map by the given scale @p factor.
  473 + ***************************************************************************/
  474 +GSkymap& GSkymap::operator*=(const double& factor)
  475 +{
  476 + // Compute total number of pixels
  477 + int n = npix() * nmaps();
  478 +
  479 + // Loop over all pixels
  480 + double* pixel = m_pixels;
  481 + for (int i = 0; i < n; ++i) {
  482 + *pixel++ *= factor;
  483 + }
  484 +
  485 + // Return this object
462 486 return *this;
463 487 }
464 488  
... ... @@ -517,7 +541,41 @@ GSkymap&amp; GSkymap::operator/=(const GSkymap&amp; map)
517 541  
518 542 } // endfor: looped over all pixels
519 543  
520   - // Return present sky map
  544 + // Return this object
  545 + return *this;
  546 +}
  547 +
  548 +
  549 +/***********************************************************************//**
  550 + * @brief Sky map division operator
  551 + *
  552 + * @param[in] factor Scale factor.
  553 + * @return Sky map.
  554 + *
  555 + * @exception GException::invalid_argument
  556 + * Division by zero error.
  557 + *
  558 + * Divides all pixels of the sky map by the given @p factor.
  559 + ***************************************************************************/
  560 +#define G_OP_UNARY_DIV2 "GSkymap::operator/=(double&)"
  561 +GSkymap& GSkymap::operator/=(const double& factor)
  562 +{
  563 + // Check for division by zero
  564 + if (factor == 0.0) {
  565 + std::string msg = "Trying to divide sky map pixels by zero.";
  566 + throw GException::invalid_argument(G_OP_UNARY_DIV2, msg);
  567 + }
  568 +
  569 + // Compute total number of pixels
  570 + int n = npix() * nmaps();
  571 +
  572 + // Loop over all pixels
  573 + double* pixel = m_pixels;
  574 + for (int i = 0; i < n; ++i) {
  575 + *pixel++ /= factor;
  576 + }
  577 +
  578 + // Return this object
521 579 return *this;
522 580 }
523 581  
... ...
test/test_GSky.cpp
  View file @961a877
... ... @@ -851,7 +851,7 @@ void TestGSky::test_GSkymap(void)
851 851 map_dst += map_src;
852 852  
853 853 // Check total in destination map. Note that the total in the destination
854   - // map should be 100 times the total in the source maps as the operator
  854 + // map should be 100 times the total in the source map as the operator
855 855 // is expected to perform strict bi-linear interpolation
856 856 double total_dst = 0.0;
857 857 for (int pix = 0; pix < map_dst.npix(); ++pix) {
... ... @@ -860,7 +860,7 @@ void TestGSky::test_GSkymap(void)
860 860 }
861 861 }
862 862 total_dst /= 100.0;
863   - test_value(total_dst, total_src, 1.0e-3, "Test operator+");
  863 + test_value(total_dst, total_src, 1.0e-3, "Test operator+(GSkymap)");
864 864  
865 865 // Subtract pixels from destination map
866 866 map_dst -= map_src;
... ... @@ -873,7 +873,55 @@ void TestGSky::test_GSkymap(void)
873 873 total_dst += map_dst(pix,k);
874 874 }
875 875 }
876   - test_value(total_dst, 0.0, 1.0e-3, "Test operator+");
  876 + test_value(total_dst, 0.0, 1.0e-3, "Test operator-(GSkymap)");
  877 +
  878 + // Check map multiplication
  879 + GSkymap test_map = map_src;
  880 + test_map *= map_src;
  881 + double total_test = 0.0;
  882 + double total_ref = 0.0;
  883 + for (int pix = 0; pix < test_map.npix(); ++pix) {
  884 + for (int k = 0; k < test_map.nmaps(); ++k) {
  885 + total_test += test_map(pix,k);
  886 + total_ref += map_src(pix,k) * map_src(pix,k);
  887 + }
  888 + }
  889 + test_value(total_test, total_ref, 1.0e-3, "Test operator*(GSkymap)");
  890 +
  891 + // Check map division
  892 + test_map = map_src;
  893 + test_map /= map_src;
  894 + total_test = 0.0;
  895 + total_ref = 0.0;
  896 + for (int pix = 0; pix < test_map.npix(); ++pix) {
  897 + for (int k = 0; k < test_map.nmaps(); ++k) {
  898 + total_test += test_map(pix,k);
  899 + total_ref += map_src(pix,k) / map_src(pix,k);
  900 + }
  901 + }
  902 + test_value(total_test, total_ref, 1.0e-3, "Test operator/(GSkymap)");
  903 +
  904 + // Check map scaling
  905 + test_map = map_src;
  906 + test_map *= 3.3;
  907 + total_test = 0.0;
  908 + for (int pix = 0; pix < test_map.npix(); ++pix) {
  909 + for (int k = 0; k < test_map.nmaps(); ++k) {
  910 + total_test += test_map(pix,k);
  911 + }
  912 + }
  913 + test_value(total_test, total_src*3.3, 1.0e-3, "Test operator*(double)");
  914 +
  915 + // Check map division
  916 + test_map = map_src;
  917 + test_map /= 3.3;
  918 + total_test = 0.0;
  919 + for (int pix = 0; pix < test_map.npix(); ++pix) {
  920 + for (int k = 0; k < test_map.nmaps(); ++k) {
  921 + total_test += test_map(pix,k);
  922 + }
  923 + }
  924 + test_value(total_test, total_src/3.3, 1.0e-3, "Test operator/(double)");
877 925  
878 926 // Save maps
879 927 map_src.save("test_map_src.fits", true);
... ...
test/test_GSky.py
  View file @961a877
... ... @@ -332,31 +332,45 @@ class Test(GPythonTestSuite):
332 332  
333 333 # Addition operator
334 334 map += map_b
335   - self.test_value(map[0], 2.0);
336   - self.test_value(map[1], 4.0);
337   - self.test_value(map[2], 6.0);
338   - self.test_value(map[3], 8.0);
  335 + self.test_value(map[0], 2.0)
  336 + self.test_value(map[1], 4.0)
  337 + self.test_value(map[2], 6.0)
  338 + self.test_value(map[3], 8.0)
339 339  
340 340 # Multiplication operator
341 341 map *= map_b
342   - self.test_value(map[0], 2.0);
343   - self.test_value(map[1], 8.0);
344   - self.test_value(map[2], 18.0);
345   - self.test_value(map[3], 32.0);
  342 + self.test_value(map[0], 2.0)
  343 + self.test_value(map[1], 8.0)
  344 + self.test_value(map[2], 18.0)
  345 + self.test_value(map[3], 32.0)
346 346  
347 347 # Subtraction operator
348 348 map -= map_b
349   - self.test_value(map[0], 1.0);
350   - self.test_value(map[1], 6.0);
351   - self.test_value(map[2], 15.0);
352   - self.test_value(map[3], 28.0);
  349 + self.test_value(map[0], 1.0)
  350 + self.test_value(map[1], 6.0)
  351 + self.test_value(map[2], 15.0)
  352 + self.test_value(map[3], 28.0)
353 353  
354 354 # Division operator
355 355 map /= map_b
356   - self.test_value(map[0], 1.0);
357   - self.test_value(map[1], 3.0);
358   - self.test_value(map[2], 5.0);
359   - self.test_value(map[3], 7.0);
  356 + self.test_value(map[0], 1.0)
  357 + self.test_value(map[1], 3.0)
  358 + self.test_value(map[2], 5.0)
  359 + self.test_value(map[3], 7.0)
  360 +
  361 + # Scaling operator
  362 + map *= 2.0
  363 + self.test_value(map[0], 2.0)
  364 + self.test_value(map[1], 6.0)
  365 + self.test_value(map[2], 10.0)
  366 + self.test_value(map[3], 14.0)
  367 +
  368 + # Division operator
  369 + map /= 2.0
  370 + self.test_value(map[0], 1.0)
  371 + self.test_value(map[1], 3.0)
  372 + self.test_value(map[2], 5.0)
  373 + self.test_value(map[3], 7.0)
360 374  
361 375 # Return
362 376 return
... ...