LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
tractor.py
Go to the documentation of this file.
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"
23 into a lsst.shapelet.MultiShapeletBasis.
24 
25 Please see the README file in the data directory of the lsst.shapelet
26 package for more information.
27 """
28 import numpy
29 import os
30 import re
31 import sys
32 import warnings
33 import pickle
34 
36 from .radialProfile import RadialProfile
37 from .multiShapeletBasis import MultiShapeletBasis
38 from .shapeletFunction import ShapeletFunction
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)
60  warnings.warn("No C++ profile for multi-Gaussian pickle file '%s'" % filename)
61  continue
62  with open(os.path.join(dataDir, filename), 'rb') as stream:
63  if sys.version_info[0] >= 3:
64  array = pickle.load(stream, encoding='latin1')
65  else:
66  array = pickle.load(stream)
67  amplitudes = array[:nComponents]
68  amplitudes /= amplitudes.sum()
69  variances = array[nComponents:]
70  if amplitudes.shape != (nComponents,) or variances.shape != (nComponents,):
71  warnings.warn("Unknown format for multi-Gaussian pickle file '%s'" % filename)
72  continue
73  basis = MultiShapeletBasis(1)
74  for amplitude, variance in zip(amplitudes, variances):
75  radius = variance**0.5
76  matrix = numpy.array([[amplitude / ShapeletFunction.FLUX_FACTOR]], dtype=float)
77  basis.addComponent(radius, 0, matrix)
78  profile.registerBasis(basis, nComponents, maxRadius)
79 
80 
81 # We register all the profiles at module import time, to allow C++ code to access all available profiles
82 # without having to later call Python code to unpickle them.
84 
85 
86 def evaluateRadial(basis, r, sbNormalize=False, doComponents=False):
87  """Plot a single-element MultiShapeletBasis as a radial profile.
88 
89  Parameters
90  ----------
91  sbNormalize : `bool`
92  `True` to normalize.
93  doComponents : `bool`
94  `True` to evaluate components.
95  """
97  coefficients = numpy.ones(1, dtype=float)
98  msf = basis.makeFunction(ellipse, coefficients)
99  ev = msf.evaluate()
100  n = 1
101  if doComponents:
102  n += len(msf.getComponents())
103  z = numpy.zeros((n,) + r.shape, dtype=float)
104  for j, x in enumerate(r):
105  z[0, j] = ev(x, 0.0)
106  if doComponents:
107  for i, sf in enumerate(msf.getComponents()):
108  evc = sf.evaluate()
109  for j, x in enumerate(r):
110  z[i+1, j] = evc(x, 0.0)
111  if sbNormalize:
112  z /= ev(1.0, 0.0)
113  return z
114 
115 
116 def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000):
117  """Integrate the profiles to compare relative fluxes between the true profiles
118  and their approximations.
119 
120  After normalizing by surface brightness at r=1 r_e, integrate the profiles to compare
121  relative fluxes between the true profiles and their approximations.
122 
123  Parameters
124  ----------
125  maxRadius : `float`, optional
126  Maximum radius to integrate the profile, in units of r_e.
127  nSteps : `int`, optional
128  Number of concrete points at which to evaluate the profile to
129  do the integration (we just use the trapezoidal rule).
130 
131  Returns
132  -------
133  fluxes : `dict` of `float` values
134  Dictionary of fluxes (``exp``, ``lux``, ``dev``, ``luv``, ``ser2``, ``ser3``,
135  ``ser5``, ``gexp``, ``glux``, ``gdev``, ``gluv``, ``gser2``, ``gser3``, ``gser5``)
136  """
137  radii = numpy.linspace(0.0, maxRadius, nSteps)
138  profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv",
139  "ser2", "ser3", "ser5")}
140  evaluated = {}
141  for name, profile in profiles.items():
142  evaluated[name] = profile.evaluate(radii)
143  basis = profile.getBasis(8)
144  evaluated["g" + name] = evaluateRadial(basis, radii, sbNormalize=True, doComponents=False)[0, :]
145  fluxes = {name: numpy.trapz(z*radii, radii) for name, z in evaluated.items()}
146  return fluxes
147 
148 
149 def plotSuite(doComponents=False):
150  """Plot all the profiles defined in this module together.
151 
152  Plot all the profiles defined in this module together: true exp and dev,
153  the SDSS softened/truncated lux and luv, and the multi-Gaussian approximations
154  to all of these.
155 
156  Parameters
157  ----------
158  doComponents : `bool`, optional
159  True, to plot the individual Gaussians that form the multi-Gaussian approximations.
160 
161  Returns
162  -------
163  figure : `matplotlib.figure.Figure`
164  Figure that contains the plot.
165  axes : `numpy.ndarray` of `matplotlib.axes.Axes`
166  A 2x4 NumPy array of matplotlib axes objects.
167  """
168  from matplotlib import pyplot
169  fig = pyplot.figure(figsize=(9, 4.7))
170  axes = numpy.zeros((2, 4), dtype=object)
171  r1 = numpy.logspace(-3, 0, 1000, base=10)
172  r2 = numpy.linspace(1, 10, 1000)
173  r = [r1, r2]
174  for i in range(2):
175  for j in range(4):
176  axes[i, j] = fig.add_subplot(2, 4, i*4+j+1)
177  profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv")}
178  basis = {name: profiles[name].getBasis(8) for name in profiles}
179  z = numpy.zeros((2, 4), dtype=object)
180  colors = ("k", "g", "b", "r")
181  fig.subplots_adjust(wspace=0.025, hspace=0.025, bottom=0.15, left=0.1, right=0.98, top=0.92)
182  centers = [None, None]
183  for i in range(2): # 0=profile, 1=relative error
184  for j in range(0, 4, 2): # grid columns: 0=exp-like, 2=dev-like
185  bbox0 = axes[i, j].get_position()
186  bbox1 = axes[i, j+1].get_position()
187  bbox1.x0 = bbox0.x1 - 0.06
188  bbox0.x1 = bbox1.x0
189  centers[j//2] = 0.5*(bbox0.x0 + bbox1.x1)
190  axes[i, j].set_position(bbox0)
191  axes[i, j+1].set_position(bbox1)
192  for j in range(0, 2):
193  z[0, j] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
194  for k in ("exp", "lux")]
195  z[0, j][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("exp", "lux")]
196  z[0, j+2] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
197  for k in ("dev", "luv")]
198  z[0, j+2][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("dev", "luv")]
199  methodNames = [["loglog", "semilogy"], ["semilogx", "plot"]]
200  for j in range(0, 4): # grid columns
201  z[1, j] = [(z[0, j][0][0, :] - z[0, j][i][0, :])/z[0, j][0][0, :] for i in range(0, 4)]
202  handles = []
203  method0 = getattr(axes[0, j], methodNames[0][j%2])
204  method1 = getattr(axes[1, j], methodNames[1][j%2])
205  for k in range(4):
206  y0 = z[0, j][k]
207  handles.append(method0(r[j%2], y0[0, :], color=colors[k])[0])
208  if doComponents:
209  for ll in range(1, y0.shape[0]):
210  method0(r[j%2], y0[ll, :], color=colors[k], alpha=0.25)
211  method1(r[j%2], z[1, j][k], color=colors[k])
212  axes[0, j].set_xticklabels([])
213  axes[0, j].set_ylim(1E-6, 1E3)
214  axes[1, j].set_ylim(-0.2, 1.0)
215  for i, label in enumerate(("profile", "relative error")):
216  axes[i, 0].set_ylabel(label)
217  for t in axes[i, 0].get_yticklabels():
218  t.set_fontsize(11)
219  for j in range(1, 4):
220  axes[0, j].set_yticklabels([])
221  axes[1, j].set_yticklabels([])
222  xticks = [['$\\mathdefault{10^{%d}}$' % i for i in range(-3, 1)],
223  [str(i) for i in range(1, 11)]]
224  xticks[0][-1] = ""
225  xticks[1][-1] = ""
226  for j in range(0, 4):
227  axes[1, j].set_xticklabels(xticks[j%2])
228  for t in axes[1, j].get_xticklabels():
229  t.set_fontsize(11)
230  fig.legend(handles, ["exp/dev", "lux/luv", "approx exp/dev", "approx lux/luv"],
231  loc='lower center', ncol=4)
232  fig.text(centers[0], 0.95, "exponential", ha='center', weight='bold')
233  fig.text(centers[1], 0.95, "de Vaucouleur", ha='center', weight='bold')
234  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
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000)
Definition: tractor.py:116
def evaluateRadial(basis, r, sbNormalize=False, doComponents=False)
Definition: tractor.py:86
def registerRadialProfiles()
Definition: tractor.py:41
def plotSuite(doComponents=False)
Definition: tractor.py:149