LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
tractor.py
Go to the documentation of this file.
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.
24
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
33
34import lsst.afw.geom
35
36from ._shapeletLib import RadialProfile, MultiShapeletBasis, ShapeletFunction
37
38_LOG = logging.getLogger(__name__)
39
40
42 """Register the pickled profiles in the data directory with the RadialProfile
43 singleton registry.
44
45 This should only be called at import time by this module; it's only a function to
46 avoid polluting the module namespace with all the local variables used here.
47 """
48 dataDir = os.path.join(os.environ["SHAPELET_DIR"], "data")
49 regex = re.compile(r"([a-z]+\d?)_K(\d+)_MR(\d+)\.pickle")
50 for filename in os.listdir(dataDir):
51 match = regex.match(filename)
52 if not match:
53 continue
54 name = match.group(1)
55 nComponents = int(match.group(2))
56 maxRadius = int(match.group(3))
57 try:
58 profile = RadialProfile.get(name)
59 except Exception:
60 _LOG.warning("No C++ profile for multi-Gaussian pickle file '%s'", filename)
61 continue
62 with open(os.path.join(dataDir, filename), 'rb') as stream:
63 array = pickle.load(stream, encoding='latin1')
64 amplitudes = array[:nComponents]
65 amplitudes /= amplitudes.sum()
66 variances = array[nComponents:]
67 if amplitudes.shape != (nComponents,) or variances.shape != (nComponents,):
68 _LOG.warning("Unknown format for multi-Gaussian pickle file '%s'", filename)
69 continue
70 basis = MultiShapeletBasis(1)
71 for amplitude, variance in zip(amplitudes, variances):
72 radius = variance**0.5
73 matrix = numpy.array([[amplitude / ShapeletFunction.FLUX_FACTOR]], dtype=float)
74 basis.addComponent(radius, 0, matrix)
75 profile.registerBasis(basis, nComponents, maxRadius)
76
77
78# We register all the profiles at module import time, to allow C++ code to access all available profiles
79# without having to later call Python code to unpickle them.
81
82
83def evaluateRadial(basis, r, sbNormalize=False, doComponents=False):
84 """Plot a single-element MultiShapeletBasis as a radial profile.
85
86 Parameters
87 ----------
88 sbNormalize : `bool`
89 `True` to normalize.
90 doComponents : `bool`
91 `True` to evaluate components.
92 """
94 coefficients = numpy.ones(1, dtype=float)
95 msf = basis.makeFunction(ellipse, coefficients)
96 ev = msf.evaluate()
97 n = 1
98 if doComponents:
99 n += len(msf.getComponents())
100 z = numpy.zeros((n,) + r.shape, dtype=float)
101 for j, x in enumerate(r):
102 z[0, j] = ev(x, 0.0)
103 if doComponents:
104 for i, sf in enumerate(msf.getComponents()):
105 evc = sf.evaluate()
106 for j, x in enumerate(r):
107 z[i+1, j] = evc(x, 0.0)
108 if sbNormalize:
109 z /= ev(1.0, 0.0)
110 return z
111
112
113def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000):
114 """Integrate the profiles to compare relative fluxes between the true profiles
115 and their approximations.
116
117 After normalizing by surface brightness at r=1 r_e, integrate the profiles to compare
118 relative fluxes between the true profiles and their approximations.
119
120 Parameters
121 ----------
122 maxRadius : `float`, optional
123 Maximum radius to integrate the profile, in units of r_e.
124 nSteps : `int`, optional
125 Number of concrete points at which to evaluate the profile to
126 do the integration (we just use the trapezoidal rule).
127
128 Returns
129 -------
130 fluxes : `dict` of `float` values
131 Dictionary of fluxes (``exp``, ``lux``, ``dev``, ``luv``, ``ser2``, ``ser3``,
132 ``ser5``, ``gexp``, ``glux``, ``gdev``, ``gluv``, ``gser2``, ``gser3``, ``gser5``)
133 """
134 radii = numpy.linspace(0.0, maxRadius, nSteps)
135 profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv",
136 "ser2", "ser3", "ser5")}
137 evaluated = {}
138 for name, profile in profiles.items():
139 evaluated[name] = profile.evaluate(radii)
140 basis = profile.getBasis(8)
141 evaluated["g" + name] = evaluateRadial(basis, radii, sbNormalize=True, doComponents=False)[0, :]
142 fluxes = {name: numpy.trapz(z*radii, radii) for name, z in evaluated.items()}
143 return fluxes
144
145
146def plotSuite(doComponents=False):
147 """Plot all the profiles defined in this module together.
148
149 Plot all the profiles defined in this module together: true exp and dev,
150 the SDSS softened/truncated lux and luv, and the multi-Gaussian approximations
151 to all of these.
152
153 Parameters
154 ----------
155 doComponents : `bool`, optional
156 True, to plot the individual Gaussians that form the multi-Gaussian approximations.
157
158 Returns
159 -------
160 figure : `matplotlib.figure.Figure`
161 Figure that contains the plot.
162 axes : `numpy.ndarray` of `matplotlib.axes.Axes`
163 A 2x4 NumPy array of matplotlib axes objects.
164 """
165 from matplotlib import pyplot
166 fig = pyplot.figure(figsize=(9, 4.7))
167 axes = numpy.zeros((2, 4), dtype=object)
168 r1 = numpy.logspace(-3, 0, 1000, base=10)
169 r2 = numpy.linspace(1, 10, 1000)
170 r = [r1, r2]
171 for i in range(2):
172 for j in range(4):
173 axes[i, j] = fig.add_subplot(2, 4, i*4+j+1)
174 profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv")}
175 basis = {name: profiles[name].getBasis(8) for name in profiles}
176 z = numpy.zeros((2, 4), dtype=object)
177 colors = ("k", "g", "b", "r")
178 fig.subplots_adjust(wspace=0.025, hspace=0.025, bottom=0.15, left=0.1, right=0.98, top=0.92)
179 centers = [None, None]
180 for i in range(2): # 0=profile, 1=relative error
181 for j in range(0, 4, 2): # grid columns: 0=exp-like, 2=dev-like
182 bbox0 = axes[i, j].get_position()
183 bbox1 = axes[i, j+1].get_position()
184 bbox1.x0 = bbox0.x1 - 0.06
185 bbox0.x1 = bbox1.x0
186 centers[j//2] = 0.5*(bbox0.x0 + bbox1.x1)
187 axes[i, j].set_position(bbox0)
188 axes[i, j+1].set_position(bbox1)
189 for j in range(0, 2):
190 z[0, j] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
191 for k in ("exp", "lux")]
192 z[0, j][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("exp", "lux")]
193 z[0, j+2] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
194 for k in ("dev", "luv")]
195 z[0, j+2][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("dev", "luv")]
196 methodNames = [["loglog", "semilogy"], ["semilogx", "plot"]]
197 for j in range(0, 4): # grid columns
198 z[1, j] = [(z[0, j][0][0, :] - z[0, j][i][0, :])/z[0, j][0][0, :] for i in range(0, 4)]
199 handles = []
200 method0 = getattr(axes[0, j], methodNames[0][j%2])
201 method1 = getattr(axes[1, j], methodNames[1][j%2])
202 for k in range(4):
203 y0 = z[0, j][k]
204 handles.append(method0(r[j%2], y0[0, :], color=colors[k])[0])
205 if doComponents:
206 for ll in range(1, y0.shape[0]):
207 method0(r[j%2], y0[ll, :], color=colors[k], alpha=0.25)
208 method1(r[j%2], z[1, j][k], color=colors[k])
209 axes[0, j].set_xticklabels([])
210 axes[0, j].set_ylim(1E-6, 1E3)
211 axes[1, j].set_ylim(-0.2, 1.0)
212 for i, label in enumerate(("profile", "relative error")):
213 axes[i, 0].set_ylabel(label)
214 for t in axes[i, 0].get_yticklabels():
215 t.set_fontsize(11)
216 for j in range(1, 4):
217 axes[0, j].set_yticklabels([])
218 axes[1, j].set_yticklabels([])
219 xticks = [['$\\mathdefault{10^{%d}}$' % i for i in range(-3, 1)],
220 [str(i) for i in range(1, 11)]]
221 xticks[0][-1] = ""
222 xticks[1][-1] = ""
223 for j in range(0, 4):
224 axes[1, j].set_xticklabels(xticks[j%2])
225 for t in axes[1, j].get_xticklabels():
226 t.set_fontsize(11)
227 fig.legend(handles, ["exp/dev", "lux/luv", "approx exp/dev", "approx lux/luv"],
228 loc='lower center', ncol=4)
229 fig.text(centers[0], 0.95, "exponential", ha='center', weight='bold')
230 fig.text(centers[1], 0.95, "de Vaucouleur", ha='center', weight='bold')
231 return fig, axes
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
Definition Axes.h:47
An ellipse defined by an arbitrary BaseCore and a center point.
Definition Ellipse.h:51
evaluateRadial(basis, r, sbNormalize=False, doComponents=False)
Definition tractor.py:83
plotSuite(doComponents=False)
Definition tractor.py:146
integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000)
Definition tractor.py:113