U/danielsf/rmjarvis/cut gsco

python/lsst/sims/GalSimInterface/galSimCelestialObject.py

@@ -1,5 +1,5 @@
 from builtins import object
-import numpy
+import numpy as np
 from lsst.sims.utils import arcsecFromRadians
 
 __all__ = ["GalSimCelestialObject"]
@@ -70,34 +70,53 @@ def __init__(self, galSimType, xPupil, yPupil,
         """
         self._uniqueId = uniqueId
         self._galSimType = galSimType
-        self._xPupilRadians = xPupil
-        self._xPupilArcsec = arcsecFromRadians(xPupil)
-        self._yPupilRadians = yPupil
-        self._yPupilArcsec = arcsecFromRadians(yPupil)
-        self._halfLightRadiusRadians = halfLightRadius
-        self._halfLightRadiusArcsec = arcsecFromRadians(halfLightRadius)
-        self._minorAxisRadians = minorAxis
-        self._majorAxisRadians = majorAxis
-        self._positionAngleRadians = positionAngle
-        self._sindex = sindex
-        self._npoints = npoints
         self._fits_image_file = fits_image_file
-        self._pixel_scale = pixel_scale
-        self._rotation_angle = rotation_angle
+
         # The galsim.lens(...) function wants to be passed reduced
         # shears and magnification, so convert the WL parameters as
         # defined in phosim instance catalogs to these values.  See
         # https://github.com/GalSim-developers/GalSim/blob/releases/1.4/doc/GalSim_Quick_Reference.pdf
         # and Hoekstra, 2013, http://lanl.arxiv.org/abs/1312.5981
-        self._g1 = gamma1/(1. - kappa)   # real part of reduced shear
-        self._g2 = gamma2/(1. - kappa)   # imaginary part of reduced shear
-        self._mu = 1./((1. - kappa)**2 - (gamma1**2 + gamma2**2)) # magnification
+        g1 = gamma1/(1. - kappa)   # real part of reduced shear
+        g2 = gamma2/(1. - kappa)   # imaginary part of reduced shear
+        mu = 1./((1. - kappa)**2 - (gamma1**2 + gamma2**2)) # magnification
 
         self._fluxDict = {}
         self._sed = sed
         self._bp_dict = bp_dict
         self._photParams = photParams
 
+        # Put all the float values into a numpy array for better memory usage.
+        # (Otherwise, python allocates a shared pointer for each one of these 15 values, which
+        # adds up to a significant memory overhead.)
+        self._float_values = np.array([
+                                        xPupil,                     # xPupilRadians, 0
+                                        arcsecFromRadians(xPupil),  # xPupilArcsec, 1
+                                        yPupil,                     # yPupilRadians, 2
+                                        arcsecFromRadians(yPupil),  # yPupilArcsec, 3
+                                        halfLightRadius,            # halfLightRadiusRadians, 4
+                                        arcsecFromRadians(halfLightRadius), # halfLightRadiusArcsec, 5
+                                        minorAxis,                  # minorAxisRadians, 6
+                                        majorAxis,                  # majorAxisRadians, 7
+                                        positionAngle,              # positionAngleRadians, 8
+                                        sindex,                     # sindex, 9
+                                        pixel_scale,                # pixel_scale, 10
+                                        rotation_angle,             # rotation_angle, 11
+                                        g1,                         # g1, 12
+                                        g2,                         # g2, 13
+                                        mu,                         # mu, 14
+                                        npoints                     # npoints, 15
+                                      ], dtype=float)
+
+        # XXX: We could probably get away with np.float32 for these, but the main
+        #      advantage is from only allocating the actual memory, and not the python
+        #      pointers to the memory.  So not so much more gain to be had from switching
+        #      to single precision.
+
+    #
+    # First properties for all the non-float values.
+    #
+
     @property
     def sed(self):
         return self._sed
@@ -106,6 +125,9 @@ def sed(self):
     def uniqueId(self):
         return self._uniqueId
 
+    # XXX: I'd recommend removing all these setters and just let python raise an
+    #      AttributeError if people try to set the values, rather than raising a
+    #      RuntimeError manually.  -- MJ
     @uniqueId.setter
     def uniqueId(self, value):
         raise RuntimeError("You should not be setting the unique id on the fly; " \
@@ -120,9 +142,32 @@ def galSimType(self, value):
         raise RuntimeError("You should not be setting galSimType on the fly; " \
                            + "just instantiate a new GalSimCelestialObject")
 
+    @property
+    def npoints(self):
+        return int(self._float_values[15]+0.5)
+
+    @npoints.setter
+    def npoints(self, value):
+        raise RuntimeError("You should not be setting npoints on the fly; " \
+        + "just instantiate a new GalSimCelestialObject")
+
+    @property
+    def fits_image_file(self):
+        return self._fits_image_file
+
+    @fits_image_file.setter
+    def fits_image_file(self, value):
+        raise RuntimeError("You should not be setting fits_image_file on the fly; "
+                           "just instantiate a new GalSimCelestialObject")
+
+
+    #
+    # The float values are accessed from the numpy array called self._float_values.
+    #
+
     @property
     def xPupilRadians(self):
-        return self._xPupilRadians
+        return self._float_values[0]
 
     @xPupilRadians.setter
     def xPupilRadians(self, value):
@@ -132,7 +177,7 @@ def xPupilRadians(self, value):
 
     @property
     def xPupilArcsec(self):
-        return self._xPupilArcsec
+        return self._float_values[1]
 
     @xPupilArcsec.setter
     def xPupilArcsec(self, value):
@@ -142,7 +187,7 @@ def xPupilArcsec(self, value):
 
     @property
     def yPupilRadians(self):
-        return self._yPupilRadians
+        return self._float_values[2]
 
     @yPupilRadians.setter
     def yPupilRadians(self, value):
@@ -152,7 +197,7 @@ def yPupilRadians(self, value):
 
     @property
     def yPupilArcsec(self):
-        return self._yPupilArcsec
+        return self._float_values[3]
 
     @yPupilArcsec.setter
     def yPupilArcsec(self, value):
@@ -162,7 +207,7 @@ def yPupilArcsec(self, value):
 
     @property
     def halfLightRadiusRadians(self):
-        return self._halfLightRadiusRadians
+        return self._float_values[4]
 
     @halfLightRadiusRadians.setter
     def halfLightRadiusRadians(self, value):
@@ -172,7 +217,7 @@ def halfLightRadiusRadians(self, value):
 
     @property
     def halfLightRadiusArcsec(self):
-        return self._halfLightRadiusArcsec
+        return self._float_values[5]
 
     @halfLightRadiusArcsec.setter
     def halfLightRadiusArcsec(self, value):
@@ -182,7 +227,7 @@ def halfLightRadiusArcsec(self, value):
 
     @property
     def minorAxisRadians(self):
-        return self._minorAxisRadians
+        return self._float_values[6]
 
     @minorAxisRadians.setter
     def minorAxisRadians(self, value):
@@ -192,7 +237,7 @@ def minorAxisRadians(self, value):
 
     @property
     def majorAxisRadians(self):
-        return self._majorAxisRadians
+        return self._float_values[7]
 
     @majorAxisRadians.setter
     def majorAxisRadians(self, value):
@@ -202,7 +247,7 @@ def majorAxisRadians(self, value):
 
     @property
     def positionAngleRadians(self):
-        return self._positionAngleRadians
+        return self._float_values[8]
 
     @positionAngleRadians.setter
     def positionAngleRadians(self, value):
@@ -212,35 +257,16 @@ def positionAngleRadians(self, value):
 
     @property
     def sindex(self):
-        return self._sindex
+        return self._float_values[9]
 
     @sindex.setter
     def sindex(self, value):
         raise RuntimeError("You should not be setting sindex on the fly; " \
         + "just instantiate a new GalSimCelestialObject")
 
-
-    @property
-    def npoints(self):
-        return self._npoints
-
-    @npoints.setter
-    def npoints(self, value):
-        raise RuntimeError("You should not be setting npoints on the fly; " \
-        + "just instantiate a new GalSimCelestialObject")
-
-    @property
-    def fits_image_file(self):
-        return self._fits_image_file
-
-    @fits_image_file.setter
-    def fits_image_file(self, value):
-        raise RuntimeError("You should not be setting fits_image_file on the fly; "
-                           "just instantiate a new GalSimCelestialObject")
-
     @property
     def pixel_scale(self):
-        return self._pixel_scale
+        return self._float_values[10]
 
     @pixel_scale.setter
     def pixel_scale(self, value):
@@ -249,7 +275,7 @@ def pixel_scale(self, value):
 
     @property
     def rotation_angle(self):
-        return self._rotation_angle
+        return self._float_values[11]
 
     @rotation_angle.setter
     def rotation_angle(self, value):
@@ -258,7 +284,7 @@ def rotation_angle(self, value):
 
     @property
     def g1(self):
-        return self._g1
+        return self._float_values[12]
 
     @g1.setter
     def g1(self, value):
@@ -268,7 +294,7 @@ def g1(self, value):
 
     @property
     def g2(self):
-        return self._g2
+        return self._float_values[13]
 
     @g2.setter
     def g2(self, value):
@@ -278,14 +304,15 @@ def g2(self, value):
 
     @property
     def mu(self):
-        return self._mu
+        return self._float_values[14]
 
     @mu.setter
     def mu(self, value):
         raise RuntimeError("You should not be setting mu on the fly; " \
         + "just instantiate a new GalSimCelestialObject")
 
 
+
     def flux(self, band):
         """
         @param [in] band is the name of a bandpass

tests/test_celestial_object.py

@@ -0,0 +1,99 @@
+import unittest
+import numpy as np
+import lsst.utils.tests
+from lsst.sims.utils import arcsecFromRadians
+from lsst.sims.GalSimInterface import GalSimCelestialObject
+from lsst.sims.photUtils import Sed, BandpassDict, Bandpass
+from lsst.sims.photUtils import PhotometricParameters
+
+
+def setup_module(module):
+    lsst.utils.tests.init()
+
+
+class GSCOTestCase(unittest.TestCase):
+
+    def test_getters(self):
+        """
+        Verify that what goes into __init__ comes out of the getters
+        """
+        xpupil_rad = 2.176
+        ypupil_rad = 3.2112
+        hlr = 1.632
+        minor_axis_rad = 9.124
+        major_axis_rad = 21.2684
+        position_angle_rad = 0.1334
+        sindex = 4.387
+        npoints = 19
+        pixel_scale=25.134
+        rotation_angle_rad = 0.9335
+        gamma1 = 0.512
+        gamma2 = 3.4456
+        kappa = 1.747
+        phot_params = PhotometricParameters()
+
+        rng = np.random.RandomState(88)
+        wav = np.arange(0.1, 200.0, 0.17)
+        spec = Sed(wavelen=wav, flambda=rng.random_sample(len(wav)))
+
+        # copy spec for comparison below
+        spec2 = Sed(wavelen=spec.wavelen, flambda=spec.flambda)
+
+        # keep BandpassDict on the same wavelength grid as Sed,
+        # otherwise, the random initialization results in flux==NaN
+        bp_list = []
+        bp_name_list = []
+        for bp_name in 'abcd':
+            sb = rng.random_sample(len(wav))
+            bp = Bandpass(wavelen=wav, sb=sb)
+            bp_list.append(bp)
+            bp_name_list.append(bp_name)
+
+        bp_dict = BandpassDict(bp_list, bp_name_list)
+
+        gso = GalSimCelestialObject('pointSource',
+                                    xpupil_rad, ypupil_rad,
+                                    hlr, minor_axis_rad, major_axis_rad,
+                                    position_angle_rad, sindex,
+                                    spec, bp_dict, phot_params,
+                                    npoints, 'bob', pixel_scale,
+                                    rotation_angle_rad,
+                                    gamma1=gamma1, gamma2=gamma2,
+                                    kappa=kappa, uniqueId=111)
+
+        self.assertAlmostEqual(gso.xPupilRadians/xpupil_rad, 1.0, 10)
+        self.assertAlmostEqual(gso.xPupilArcsec/arcsecFromRadians(xpupil_rad), 1.0, 10)
+        self.assertAlmostEqual(gso.yPupilRadians/ypupil_rad, 1.0, 10)
+        self.assertAlmostEqual(gso.yPupilArcsec/arcsecFromRadians(ypupil_rad), 1.0, 10)
+        self.assertAlmostEqual(gso.halfLightRadiusRadians/hlr, 1.0, 10)
+        self.assertAlmostEqual(gso.halfLightRadiusArcsec/arcsecFromRadians(hlr), 1.0, 10)
+        self.assertEqual(gso.uniqueId, 111)
+        self.assertEqual(gso.galSimType, 'pointSource')
+        self.assertEqual(gso.npoints, npoints)
+        self.assertAlmostEqual(gso.minorAxisRadians/minor_axis_rad, 1.0, 10)
+        self.assertAlmostEqual(gso.majorAxisRadians/major_axis_rad, 1.0, 10)
+        self.assertAlmostEqual(gso.positionAngleRadians/position_angle_rad, 1.0, 10)
+        self.assertAlmostEqual(gso.sindex/sindex, 1.0, 10)
+        self.assertAlmostEqual(gso.pixel_scale/pixel_scale, 1.0, 10)
+        self.assertAlmostEqual(gso.rotation_angle/rotation_angle_rad, 1.0, 10)
+        g1 = gamma1/(1.0-kappa)
+        self.assertAlmostEqual(gso.g1/g1, 1.0, 10)
+        g2 = gamma2/(1.0-kappa)
+        self.assertAlmostEqual(gso.g2/g2, 1.0, 10)
+        mu = 1.0/((1.0-kappa)**2-(gamma1**2+gamma2**2))
+        self.assertAlmostEqual(gso.mu/mu, 1.0, 10)
+        self.assertEqual(gso.fits_image_file, 'bob')
+        self.assertEqual(gso.sed, spec2)
+        for bp_name in bp_dict:
+            ff = spec2.calcADU(bp_dict[bp_name], phot_params)
+            ff *= phot_params.gain
+            self.assertTrue(np.isfinite(ff))
+            self.assertAlmostEqual(gso.flux(bp_name)/ff, 1.0, 10)
+
+
+class MemoryTestClass(lsst.utils.tests.MemoryTestCase):
+    pass
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()