DM-7256: Python 3 for Shapelet

examples/shapeletBases.py

@@ -21,20 +21,24 @@
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
 
+from __future__ import print_function
+from builtins import range
 import lsst.shapelet
 from lsst.afw import geom
 from lsst.afw.geom import ellipses
 from matplotlib import pyplot
 from matplotlib import ticker
 import numpy
 
+
 def makeBasisImages(basis, x, y):
     z = numpy.zeros((y.size, x.size, lsst.shapelet.computeSize(basis.getOrder())), dtype=float)
     for i, py in enumerate(y):
         for j, px in enumerate(x):
-            basis.fillEvaluation(z[i,j,:], float(px), float(py))
+            basis.fillEvaluation(z[i, j, :], float(px), float(py))
     return z
 
+
 def compareMoments(basis, x, y, z):
     e = basis.evaluate()
     monopole = z.sum()
@@ -45,9 +49,10 @@ def compareMoments(basis, x, y, z):
         (dx**2 * z).sum() / monopole,
         (dy**1 * z).sum() / monopole,
         (dx * dy * z).sum() / monopole
-        )
-    print ellipses.Ellipse(quadrupole, monopole)
-    print e.computeMoments()
+    )
+    print(ellipses.Ellipse(quadrupole, monopole))
+    print(e.computeMoments())
+
 
 def plotBasisImages(basis, z):
     n = basis.getOrder()
@@ -58,7 +63,7 @@ def plotBasisImages(basis, z):
     for i in range(n+1):
         for j in range(i+1):
             axes = pyplot.subplot(n+1, n+1, (n+1) * i + j + 1)
-            axes.imshow(z[:,:,k], interpolation='nearest', origin='lower', vmin=vmin, vmax=vmax)
+            axes.imshow(z[:, :, k], interpolation='nearest', origin='lower', vmin=vmin, vmax=vmax)
             axes.yaxis.set_major_locator(ticker.NullLocator())
             axes.xaxis.set_major_locator(ticker.NullLocator())
             if basis.getBasisType() == lsst.shapelet.HERMITE:
@@ -67,10 +72,11 @@ def plotBasisImages(basis, z):
                 pyplot.xlabel("p=%d, q=%d (%s)" % (i-j/2, j/2, "Im" if j % 2 else "Re"))
             k += 1
 
+
 def processBasis(basis, x, y):
     z = makeBasisImages(basis, x, y)
     plotBasisImages(basis, z)
-    
+
 
 def main():
     x = numpy.linspace(-5, 5, 101)

examples/shapeletConvolve.py

@@ -28,16 +28,18 @@
 from matplotlib import ticker
 import numpy
 
+
 def plotShapeletFunction(axes, func, x, y):
     z = numpy.zeros((y.size, x.size), dtype=float)
     ev = func.evaluate()
     for i, py in enumerate(y):
         for j, px in enumerate(x):
-            z[i,j] = ev(float(px), float(py))
+            z[i, j] = ev(float(px), float(py))
     axes.imshow(z, interpolation='nearest', origin='lower')
     axes.yaxis.set_major_locator(ticker.NullLocator())
     axes.xaxis.set_major_locator(ticker.NullLocator())
 
+
 def main():
     x = numpy.linspace(-5, 5, 101)
     y = numpy.linspace(-5, 5, 101)

python/lsst/__init__.py

@@ -1,2 +1,3 @@
-import pkgutil, lsstimport
+import pkgutil
+import lsstimport
 __path__ = pkgutil.extend_path(__path__, __name__)

python/lsst/shapelet/__init__.py

@@ -1,29 +1,33 @@
-# 
+from __future__ import absolute_import
+from builtins import range
+from builtins import object
+#
 # LSST Data Management System
 # Copyright 2008, 2009, 2010 LSST Corporation.
-# 
+#
 # This product includes software developed by the
 # LSST Project (http://www.lsst.org/).
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
-# 
+#
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
-# 
-# You should have received a copy of the LSST License Statement and 
-# the GNU General Public License along with this program.  If not, 
+#
+# You should have received a copy of the LSST License Statement and
+# the GNU General Public License along with this program.  If not,
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
 
-from shapeletLib import *
+from .shapeletLib import *
 from . import tractor
 from lsst.afw.geom.ellipses import Quadrupole as EllipseCore
 
+
 class IndexGenerator(object):
     """
     Base class for shapelet index generators.
@@ -45,32 +49,34 @@ def __init__(self, order):
     def __len__(self):
         return self.size
 
+
 class HermiteIndexGenerator(IndexGenerator):
     """
     Iterable that generates tuples of (i, nx, ny) in which:
      - 'i' is the overall coefficient index for a 2-d shapelet expansion (just counts from zero)
      - 'nx' is the order of the x expansion
      - 'ny' is the order of the y expansion
      """
-    
+
     def __iter__(self):
         i = 0
-        for n in xrange(0, self.order+1):
-            for nx in xrange(0, n+1):
+        for n in range(0, self.order+1):
+            for nx in range(0, n+1):
                 yield (i, nx, n - nx)
                 i += 1
 
+
 class LaguerreIndexGenerator(IndexGenerator):
     """
     Iterable that generates tuples of (i, p, q, re) in which:
      - 'i' is the overall coefficient index for a 2-d shapelet expansion (just counts from zero)
      - 'p' and 'q' are the indices of the polar shapelet expansion (see BasisTypeEnum).
      - 're' is True if this the real part of the coefficient
      """
-    
+
     def __iter__(self):
         i = 0
-        for n in xrange(0, self.order+1):
+        for n in range(0, self.order+1):
             p = n
             q = 0
             while p > q:
@@ -83,4 +89,3 @@ def __iter__(self):
             if p == q:
                 yield (i, p, q, True)
                 i += 1
-

python/lsst/shapelet/tests.py

@@ -24,6 +24,9 @@
 Test utility code for shapelets library; here so it can be used
 in multiple test scripts and tests in downstream packages.
 """
+from __future__ import print_function
+from builtins import zip
+from builtins import range
 
 import numpy
 try:
@@ -35,6 +38,7 @@
 import lsst.afw.geom
 import lsst.afw.geom.ellipses
 
+
 class ShapeletTestCase(lsst.utils.tests.TestCase):
 
     @staticmethod
@@ -49,7 +53,7 @@ def makeImage(function, x, y):
         e = function.evaluate()
         for i, py in enumerate(y):
             for j, px in enumerate(x):
-                z[i,j] = e(float(px), float(py))
+                z[i, j] = e(float(px), float(py))
         return z
 
     @staticmethod
@@ -63,9 +67,9 @@ def makeRandomShapeletFunction(order=2, zeroCenter=False, ellipse=None, scale=1.
                     float(numpy.random.uniform(low=1, high=2)),
                     float(numpy.random.uniform(low=1, high=2)),
                     float(numpy.random.uniform(low=0, high=numpy.pi))
-                    ),
+                ),
                 center
-                )
+            )
         coefficients = numpy.random.randn(lsst.shapelet.computeSize(order))
         result = lsst.shapelet.ShapeletFunction(order, lsst.shapelet.HERMITE, coefficients)
         result.setEllipse(ellipse)
@@ -119,7 +123,7 @@ def checkMoments(self, function, x, y, z):
             (gx**2 * z).sum() / m,
             (gy**2 * z).sum() / m,
             (gx * gy * z).sum() / m
-            )
+        )
         imageMoments = lsst.afw.geom.ellipses.Ellipse(quadrupole, dipole)
         shapeletMoments = function.evaluate().computeMoments()
         self.assertClose(imageMoments.getCenter().getX(), shapeletMoments.getCenter().getX(), rtol=1E-3)
@@ -147,7 +151,7 @@ def checkConvolution(self, f1, f2):
         self.assertClose(ic1.getArray(), ic2.getArray())
         out = lsst.afw.image.ImageD(bbox)
         if scipy is None:
-            print "Skipping convolution test; scipy could not be imported."
+            print("Skipping convolution test; scipy could not be imported.")
             return
         # I'm using scipy.ndimage to convolve test images, because I can't figure
         # out how to make afw do it (afw can convolve images with kernels, but two similarly-sized

python/lsst/shapelet/tractor.py

@@ -26,18 +26,25 @@
 Please see the README file in the data directory of the lsst.shapelet
 package for more information.
 """
+from future import standard_library
+standard_library.install_aliases()
+from builtins import zip
+from builtins import str
+from builtins import range
 
 import numpy
 import os
 import re
+import sys
 import warnings
 try:
-    import cPickle as pickle
+    import pickle as pickle
 except ImportError:
     import pickle
 
 from .shapeletLib import *
 
+
 def registerRadialProfiles():
     """Register the pickled profiles in the data directory with the RadialProfile singleton registry.
 
@@ -48,7 +55,8 @@ def registerRadialProfiles():
     regex = re.compile(r"([a-z]+\d?)_K(\d+)_MR(\d+)\.pickle")
     for filename in os.listdir(dataDir):
         match = regex.match(filename)
-        if not match: continue
+        if not match:
+            continue
         name = match.group(1)
         nComponents = int(match.group(2))
         maxRadius = int(match.group(3))
@@ -57,8 +65,11 @@ def registerRadialProfiles():
         except lsst.pex.exceptions.Exception:
             warnings.warn("No C++ profile for multi-Gaussian pickle file '%s'" % filename)
             continue
-        with open(os.path.join(dataDir, filename), 'r') as stream:
-            array = pickle.load(stream)
+        with open(os.path.join(dataDir, filename), 'rb') as stream:
+            if sys.version_info[0] >= 3:
+                array = pickle.load(stream, encoding='latin1')
+            else:
+                array = pickle.load(stream)
         amplitudes = array[:nComponents]
         amplitudes /= amplitudes.sum()
         variances = array[nComponents:]
@@ -75,6 +86,7 @@ def registerRadialProfiles():
 # without having to later call Python code to unpickle them.
 registerRadialProfiles()
 
+
 def evaluateRadial(basis, r, sbNormalize=False, doComponents=False):
     """Plot a single-element MultiShapeletBasis as a radial profile.
     """
@@ -87,16 +99,17 @@ def evaluateRadial(basis, r, sbNormalize=False, doComponents=False):
         n += len(msf.getComponents())
     z = numpy.zeros((n,) + r.shape, dtype=float)
     for j, x in enumerate(r):
-        z[0,j] = ev(x, 0.0)
+        z[0, j] = ev(x, 0.0)
     if doComponents:
         for i, sf in enumerate(msf.getComponents()):
             evc = sf.evaluate()
             for j, x in enumerate(r):
-                z[i+1,j] = evc(x, 0.0)
+                z[i+1, j] = evc(x, 0.0)
     if sbNormalize:
         z /= ev(1.0, 0.0)
     return z
 
+
 def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000):
     """!
     After normalizing by surface brightness at r=1 r_e, integrate the profiles to compare
@@ -110,11 +123,11 @@ def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000):
     profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv",
                                                            "ser2", "ser3", "ser5")}
     evaluated = {}
-    for name, profile in profiles.iteritems():
+    for name, profile in profiles.items():
         evaluated[name] = profile.evaluate(radii)
         basis = profile.getBasis(8)
-        evaluated["g" + name] = evaluateRadial(basis, radii, sbNormalize=True, doComponents=False)[0,:]
-    fluxes = {name: numpy.trapz(z*radii, radii) for name, z in evaluated.iteritems()}
+        evaluated["g" + name] = evaluateRadial(basis, radii, sbNormalize=True, doComponents=False)[0, :]
+    fluxes = {name: numpy.trapz(z*radii, radii) for name, z in evaluated.items()}
     return fluxes
 
 
@@ -135,10 +148,10 @@ def plotSuite(doComponents=False):
     r = [r1, r2]
     for i in range(2):
         for j in range(4):
-            axes[i,j] = fig.add_subplot(2, 4, i*4+j+1)
+            axes[i, j] = fig.add_subplot(2, 4, i*4+j+1)
     profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv")}
     basis = {name: profiles[name].getBasis(8) for name in profiles}
-    z = numpy.zeros((2,4), dtype=object)
+    z = numpy.zeros((2, 4), dtype=object)
     colors = ("k", "g", "b", "r")
     fig.subplots_adjust(wspace=0.025, hspace=0.025, bottom=0.15, left=0.1, right=0.98, top=0.92)
     centers = [None, None]
@@ -149,17 +162,17 @@ def plotSuite(doComponents=False):
             bbox1.x0 = bbox0.x1 - 0.06
             bbox0.x1 = bbox1.x0
             centers[j/2] = 0.5*(bbox0.x0 + bbox1.x1)
-            axes[i,j].set_position(bbox0)
-            axes[i,j+1].set_position(bbox1)
-    for j in range(0,2):
-        z[0,j] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
-                  for k in ("exp", "lux")]
-        z[0,j][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis,:] for k in ("exp", "lux")]
-        z[0,j+2] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
-                    for k in ("dev", "luv")]
-        z[0,j+2][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis,:] for k in ("dev", "luv")]
+            axes[i, j].set_position(bbox0)
+            axes[i, j+1].set_position(bbox1)
+    for j in range(0, 2):
+        z[0, j] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
+                   for k in ("exp", "lux")]
+        z[0, j][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("exp", "lux")]
+        z[0, j+2] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
+                     for k in ("dev", "luv")]
+        z[0, j+2][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("dev", "luv")]
     methodNames = [["loglog", "semilogy"], ["semilogx", "plot"]]
-    for j in range(0, 4): # grid columns
+    for j in range(0, 4):  # grid columns
         y[1, j] = [(y[0, j][0][0, :] - y[0, j][i][0, :])/y[0, j][0][0, :] for i in range(0, 4)]
         handles = []
         method0 = getattr(axes[0, j], methodNames[0][j%2])
@@ -182,7 +195,7 @@ def plotSuite(doComponents=False):
         axes[0, j].set_yticklabels([])
         axes[1, j].set_yticklabels([])
     xticks = [['$\\mathdefault{10^{%d}}$' % i for i in range(-3, 1)],
-               [str(i) for i in range(1, 11)]]
+              [str(i) for i in range(1, 11)]]
     xticks[0][-1] = ""
     xticks[1][-1] = ""
     for j in range(0, 4):