Coverage for python/lsst/shapelet/tractor.py: 34%
132 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-07 23:04 +0000
1#
2# LSST Data Management System
3# Copyright 2008-2013 LSST Corporation.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <http://www.lsstcorp.org/LegalNotices/>.
21#
22"""Code to load multi-Gaussian approximations to profiles from "The Tractor"
23into a lsst.shapelet.MultiShapeletBasis.
25Please see the README file in the data directory of the lsst.shapelet
26package for more information.
27"""
28import numpy
29import os
30import re
31import pickle
32import logging
33import importlib.resources
35import lsst.afw.geom
37from ._shapeletLib import RadialProfile, MultiShapeletBasis, ShapeletFunction
39_LOG = logging.getLogger(__name__)
42def registerRadialProfiles():
43 """Register the pickled profiles in the data directory with the RadialProfile
44 singleton registry.
46 This should only be called at import time by this module; it's only a function to
47 avoid polluting the module namespace with all the local variables used here.
48 """
49 regex = re.compile(r"([a-z]+\d?)_K(\d+)_MR(\d+)\.pickle")
50 # Locate the pickled profiles that ship inside the package. Using
51 # importlib.resources rather than the EUPS-provided $SHAPELET_DIR means
52 # this works in an eups-less (pip/conda) installation as well.
53 with importlib.resources.as_file(
54 importlib.resources.files("lsst.shapelet").joinpath("resources", "data")
55 ) as dataDir:
56 for filename in os.listdir(dataDir):
57 match = regex.match(filename)
58 if not match:
59 continue
60 name = match.group(1)
61 nComponents = int(match.group(2))
62 maxRadius = int(match.group(3))
63 try:
64 profile = RadialProfile.get(name)
65 except Exception:
66 _LOG.warning("No C++ profile for multi-Gaussian pickle file '%s'", filename)
67 continue
68 with open(os.path.join(dataDir, filename), 'rb') as stream:
69 array = pickle.load(stream, encoding='latin1')
70 amplitudes = array[:nComponents]
71 amplitudes /= amplitudes.sum()
72 variances = array[nComponents:]
73 if amplitudes.shape != (nComponents,) or variances.shape != (nComponents,): 73 ↛ 74line 73 didn't jump to line 74 because the condition on line 73 was never true
74 _LOG.warning("Unknown format for multi-Gaussian pickle file '%s'", filename)
75 continue
76 basis = MultiShapeletBasis(1)
77 for amplitude, variance in zip(amplitudes, variances):
78 radius = variance**0.5
79 matrix = numpy.array([[amplitude / ShapeletFunction.FLUX_FACTOR]], dtype=float)
80 basis.addComponent(radius, 0, matrix)
81 profile.registerBasis(basis, nComponents, maxRadius)
84# We register all the profiles at module import time, to allow C++ code to access all available profiles
85# without having to later call Python code to unpickle them.
86registerRadialProfiles()
89def evaluateRadial(basis, r, sbNormalize=False, doComponents=False):
90 """Plot a single-element MultiShapeletBasis as a radial profile.
92 Parameters
93 ----------
94 sbNormalize : `bool`
95 `True` to normalize.
96 doComponents : `bool`
97 `True` to evaluate components.
98 """
99 ellipse = lsst.afw.geom.ellipses.Ellipse(lsst.afw.geom.ellipses.Axes())
100 coefficients = numpy.ones(1, dtype=float)
101 msf = basis.makeFunction(ellipse, coefficients)
102 ev = msf.evaluate()
103 n = 1
104 if doComponents: 104 ↛ 105line 104 didn't jump to line 105 because the condition on line 104 was never true
105 n += len(msf.getComponents())
106 z = numpy.zeros((n,) + r.shape, dtype=float)
107 for j, x in enumerate(r):
108 z[0, j] = ev(x, 0.0)
109 if doComponents: 109 ↛ 110line 109 didn't jump to line 110 because the condition on line 109 was never true
110 for i, sf in enumerate(msf.getComponents()):
111 evc = sf.evaluate()
112 for j, x in enumerate(r):
113 z[i+1, j] = evc(x, 0.0)
114 if sbNormalize: 114 ↛ 116line 114 didn't jump to line 116 because the condition on line 114 was always true
115 z /= ev(1.0, 0.0)
116 return z
119def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000):
120 """Integrate the profiles to compare relative fluxes between the true profiles
121 and their approximations.
123 After normalizing by surface brightness at r=1 r_e, integrate the profiles to compare
124 relative fluxes between the true profiles and their approximations.
126 Parameters
127 ----------
128 maxRadius : `float`, optional
129 Maximum radius to integrate the profile, in units of r_e.
130 nSteps : `int`, optional
131 Number of concrete points at which to evaluate the profile to
132 do the integration (we just use the trapezoidal rule).
134 Returns
135 -------
136 fluxes : `dict` of `float` values
137 Dictionary of fluxes (``exp``, ``lux``, ``dev``, ``luv``, ``ser2``, ``ser3``,
138 ``ser5``, ``gexp``, ``glux``, ``gdev``, ``gluv``, ``gser2``, ``gser3``, ``gser5``)
139 """
140 radii = numpy.linspace(0.0, maxRadius, nSteps)
141 profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv",
142 "ser2", "ser3", "ser5")}
143 evaluated = {}
144 for name, profile in profiles.items():
145 evaluated[name] = profile.evaluate(radii)
146 basis = profile.getBasis(8)
147 evaluated["g" + name] = evaluateRadial(basis, radii, sbNormalize=True, doComponents=False)[0, :]
148 fluxes = {name: numpy.trapezoid(z*radii, radii) for name, z in evaluated.items()}
149 return fluxes
152def plotSuite(doComponents=False):
153 """Plot all the profiles defined in this module together.
155 Plot all the profiles defined in this module together: true exp and dev,
156 the SDSS softened/truncated lux and luv, and the multi-Gaussian approximations
157 to all of these.
159 Parameters
160 ----------
161 doComponents : `bool`, optional
162 True, to plot the individual Gaussians that form the multi-Gaussian approximations.
164 Returns
165 -------
166 figure : `matplotlib.figure.Figure`
167 Figure that contains the plot.
168 axes : `numpy.ndarray` of `matplotlib.axes.Axes`
169 A 2x4 NumPy array of matplotlib axes objects.
170 """
171 from matplotlib import pyplot
172 fig = pyplot.figure(figsize=(9, 4.7))
173 axes = numpy.zeros((2, 4), dtype=object)
174 r1 = numpy.logspace(-3, 0, 1000, base=10)
175 r2 = numpy.linspace(1, 10, 1000)
176 r = [r1, r2]
177 for i in range(2):
178 for j in range(4):
179 axes[i, j] = fig.add_subplot(2, 4, i*4+j+1)
180 profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv")}
181 basis = {name: profiles[name].getBasis(8) for name in profiles}
182 z = numpy.zeros((2, 4), dtype=object)
183 colors = ("k", "g", "b", "r")
184 fig.subplots_adjust(wspace=0.025, hspace=0.025, bottom=0.15, left=0.1, right=0.98, top=0.92)
185 centers = [None, None]
186 for i in range(2): # 0=profile, 1=relative error
187 for j in range(0, 4, 2): # grid columns: 0=exp-like, 2=dev-like
188 bbox0 = axes[i, j].get_position()
189 bbox1 = axes[i, j+1].get_position()
190 bbox1.x0 = bbox0.x1 - 0.06
191 bbox0.x1 = bbox1.x0
192 centers[j//2] = 0.5*(bbox0.x0 + bbox1.x1)
193 axes[i, j].set_position(bbox0)
194 axes[i, j+1].set_position(bbox1)
195 for j in range(0, 2):
196 z[0, j] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
197 for k in ("exp", "lux")]
198 z[0, j][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("exp", "lux")]
199 z[0, j+2] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
200 for k in ("dev", "luv")]
201 z[0, j+2][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("dev", "luv")]
202 methodNames = [["loglog", "semilogy"], ["semilogx", "plot"]]
203 for j in range(0, 4): # grid columns
204 z[1, j] = [(z[0, j][0][0, :] - z[0, j][i][0, :])/z[0, j][0][0, :] for i in range(0, 4)]
205 handles = []
206 method0 = getattr(axes[0, j], methodNames[0][j%2])
207 method1 = getattr(axes[1, j], methodNames[1][j%2])
208 for k in range(4):
209 y0 = z[0, j][k]
210 handles.append(method0(r[j%2], y0[0, :], color=colors[k])[0])
211 if doComponents:
212 for ll in range(1, y0.shape[0]):
213 method0(r[j%2], y0[ll, :], color=colors[k], alpha=0.25)
214 method1(r[j%2], z[1, j][k], color=colors[k])
215 axes[0, j].set_xticklabels([])
216 axes[0, j].set_ylim(1E-6, 1E3)
217 axes[1, j].set_ylim(-0.2, 1.0)
218 for i, label in enumerate(("profile", "relative error")):
219 axes[i, 0].set_ylabel(label)
220 for t in axes[i, 0].get_yticklabels():
221 t.set_fontsize(11)
222 for j in range(1, 4):
223 axes[0, j].set_yticklabels([])
224 axes[1, j].set_yticklabels([])
225 xticks = [['$\\mathdefault{10^{%d}}$' % i for i in range(-3, 1)],
226 [str(i) for i in range(1, 11)]]
227 xticks[0][-1] = ""
228 xticks[1][-1] = ""
229 for j in range(0, 4):
230 axes[1, j].set_xticklabels(xticks[j%2])
231 for t in axes[1, j].get_xticklabels():
232 t.set_fontsize(11)
233 fig.legend(handles, ["exp/dev", "lux/luv", "approx exp/dev", "approx lux/luv"],
234 loc='lower center', ncol=4)
235 fig.text(centers[0], 0.95, "exponential", ha='center', weight='bold')
236 fig.text(centers[1], 0.95, "de Vaucouleur", ha='center', weight='bold')
237 return fig, axes