LSSTApplications  18.1.0
LSSTDataManagementBasePackage
Functions
lsst.shapelet.tractor Namespace Reference

Functions

def registerRadialProfiles ()
 
def evaluateRadial (basis, r, sbNormalize=False, doComponents=False)
 
def integrateNormalizedFluxes (maxRadius=20.0, nSteps=5000)
 
def plotSuite (doComponents=False)
 

Function Documentation

◆ evaluateRadial()

def lsst.shapelet.tractor.evaluateRadial (   basis,
  r,
  sbNormalize = False,
  doComponents = False 
)
Plot a single-element MultiShapeletBasis as a radial profile.

Parameters
----------
sbNormalize : `bool`
    `True` to normalize.
doComponents : `bool`
    `True` to evaluate components.

Definition at line 92 of file tractor.py.

92 def evaluateRadial(basis, r, sbNormalize=False, doComponents=False):
93  """Plot a single-element MultiShapeletBasis as a radial profile.
94 
95  Parameters
96  ----------
97  sbNormalize : `bool`
98  `True` to normalize.
99  doComponents : `bool`
100  `True` to evaluate components.
101  """
103  coefficients = numpy.ones(1, dtype=float)
104  msf = basis.makeFunction(ellipse, coefficients)
105  ev = msf.evaluate()
106  n = 1
107  if doComponents:
108  n += len(msf.getComponents())
109  z = numpy.zeros((n,) + r.shape, dtype=float)
110  for j, x in enumerate(r):
111  z[0, j] = ev(x, 0.0)
112  if doComponents:
113  for i, sf in enumerate(msf.getComponents()):
114  evc = sf.evaluate()
115  for j, x in enumerate(r):
116  z[i+1, j] = evc(x, 0.0)
117  if sbNormalize:
118  z /= ev(1.0, 0.0)
119  return z
120 
121 
def evaluateRadial(basis, r, sbNormalize=False, doComponents=False)
Definition: tractor.py:92
An ellipse defined by an arbitrary BaseCore and a center point.
Definition: Ellipse.h:51
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a...
Definition: Axes.h:47

◆ integrateNormalizedFluxes()

def lsst.shapelet.tractor.integrateNormalizedFluxes (   maxRadius = 20.0,
  nSteps = 5000 
)
Integrate the profiles to compare relative fluxes between the true profiles
and their approximations.

After normalizing by surface brightness at r=1 r_e, integrate the profiles to compare
relative fluxes between the true profiles and their approximations.

Parameters
----------
maxRadius : `float`, optional
    Maximum radius to integrate the profile, in units of r_e.
nSteps : `int`, optional
    Number of concrete points at which to evaluate the profile to
    do the integration (we just use the trapezoidal rule).

Returns
-------
fluxes : `dict` of `float` values
    Dictionary of fluxes (``exp``, ``lux``, ``dev``, ``luv``, ``ser2``, ``ser3``,
    ``ser5``, ``gexp``,  ``glux``, ``gdev``, ``gluv``, ``gser2``, ``gser3``, ``gser5``)

Definition at line 122 of file tractor.py.

122 def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000):
123  """Integrate the profiles to compare relative fluxes between the true profiles
124  and their approximations.
125 
126  After normalizing by surface brightness at r=1 r_e, integrate the profiles to compare
127  relative fluxes between the true profiles and their approximations.
128 
129  Parameters
130  ----------
131  maxRadius : `float`, optional
132  Maximum radius to integrate the profile, in units of r_e.
133  nSteps : `int`, optional
134  Number of concrete points at which to evaluate the profile to
135  do the integration (we just use the trapezoidal rule).
136 
137  Returns
138  -------
139  fluxes : `dict` of `float` values
140  Dictionary of fluxes (``exp``, ``lux``, ``dev``, ``luv``, ``ser2``, ``ser3``,
141  ``ser5``, ``gexp``, ``glux``, ``gdev``, ``gluv``, ``gser2``, ``gser3``, ``gser5``)
142  """
143  radii = numpy.linspace(0.0, maxRadius, nSteps)
144  profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv",
145  "ser2", "ser3", "ser5")}
146  evaluated = {}
147  for name, profile in profiles.items():
148  evaluated[name] = profile.evaluate(radii)
149  basis = profile.getBasis(8)
150  evaluated["g" + name] = evaluateRadial(basis, radii, sbNormalize=True, doComponents=False)[0, :]
151  fluxes = {name: numpy.trapz(z*radii, radii) for name, z in evaluated.items()}
152  return fluxes
153 
154 
def evaluateRadial(basis, r, sbNormalize=False, doComponents=False)
Definition: tractor.py:92
def integrateNormalizedFluxes(maxRadius=20.0, nSteps=5000)
Definition: tractor.py:122

◆ plotSuite()

def lsst.shapelet.tractor.plotSuite (   doComponents = False)
Plot all the profiles defined in this module together.

Plot all the profiles defined in this module together: true exp and dev,
the SDSS softened/truncated lux and luv, and the multi-Gaussian approximations
to all of these.

Parameters
----------
doComponents : `bool`, optional
    True, to plot the individual Gaussians that form the multi-Gaussian approximations.

Returns
-------
figure : `matplotlib.figure.Figure`
    Figure that contains the plot.
axes : `numpy.ndarray` of `matplotlib.axes.Axes`
    A 2x4 NumPy array of matplotlib axes objects.

Definition at line 155 of file tractor.py.

155 def plotSuite(doComponents=False):
156  """Plot all the profiles defined in this module together.
157 
158  Plot all the profiles defined in this module together: true exp and dev,
159  the SDSS softened/truncated lux and luv, and the multi-Gaussian approximations
160  to all of these.
161 
162  Parameters
163  ----------
164  doComponents : `bool`, optional
165  True, to plot the individual Gaussians that form the multi-Gaussian approximations.
166 
167  Returns
168  -------
169  figure : `matplotlib.figure.Figure`
170  Figure that contains the plot.
171  axes : `numpy.ndarray` of `matplotlib.axes.Axes`
172  A 2x4 NumPy array of matplotlib axes objects.
173  """
174  from matplotlib import pyplot
175  fig = pyplot.figure(figsize=(9, 4.7))
176  axes = numpy.zeros((2, 4), dtype=object)
177  r1 = numpy.logspace(-3, 0, 1000, base=10)
178  r2 = numpy.linspace(1, 10, 1000)
179  r = [r1, r2]
180  for i in range(2):
181  for j in range(4):
182  axes[i, j] = fig.add_subplot(2, 4, i*4+j+1)
183  profiles = {name: RadialProfile.get(name) for name in ("exp", "lux", "dev", "luv")}
184  basis = {name: profiles[name].getBasis(8) for name in profiles}
185  z = numpy.zeros((2, 4), dtype=object)
186  colors = ("k", "g", "b", "r")
187  fig.subplots_adjust(wspace=0.025, hspace=0.025, bottom=0.15, left=0.1, right=0.98, top=0.92)
188  centers = [None, None]
189  for i in range(2): # 0=profile, 1=relative error
190  for j in range(0, 4, 2): # grid columns: 0=exp-like, 2=dev-like
191  bbox0 = axes[i, j].get_position()
192  bbox1 = axes[i, j+1].get_position()
193  bbox1.x0 = bbox0.x1 - 0.06
194  bbox0.x1 = bbox1.x0
195  centers[j/2] = 0.5*(bbox0.x0 + bbox1.x1)
196  axes[i, j].set_position(bbox0)
197  axes[i, j+1].set_position(bbox1)
198  for j in range(0, 2):
199  z[0, j] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
200  for k in ("exp", "lux")]
201  z[0, j][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("exp", "lux")]
202  z[0, j+2] = [evaluateRadial(basis[k], r[j], sbNormalize=True, doComponents=doComponents)
203  for k in ("dev", "luv")]
204  z[0, j+2][0:0] = [profiles[k].evaluate(r[j])[numpy.newaxis, :] for k in ("dev", "luv")]
205  methodNames = [["loglog", "semilogy"], ["semilogx", "plot"]]
206  for j in range(0, 4): # grid columns
207  z[1, j] = [(z[0, j][0][0, :] - z[0, j][i][0, :])/z[0, j][0][0, :] for i in range(0, 4)]
208  handles = []
209  method0 = getattr(axes[0, j], methodNames[0][j%2])
210  method1 = getattr(axes[1, j], methodNames[1][j%2])
211  for k in range(4):
212  y0 = z[0, j][k]
213  handles.append(method0(r[j%2], y0[0, :], color=colors[k])[0])
214  if doComponents:
215  for l in range(1, y0.shape[0]):
216  method0(r[j%2], y0[l, :], color=colors[k], alpha=0.25)
217  method1(r[j%2], z[1, j][k], color=colors[k])
218  axes[0, j].set_xticklabels([])
219  axes[0, j].set_ylim(1E-6, 1E3)
220  axes[1, j].set_ylim(-0.2, 1.0)
221  for i, label in enumerate(("profile", "relative error")):
222  axes[i, 0].set_ylabel(label)
223  for t in axes[i, 0].get_yticklabels():
224  t.set_fontsize(11)
225  for j in range(1, 4):
226  axes[0, j].set_yticklabels([])
227  axes[1, j].set_yticklabels([])
228  xticks = [['$\\mathdefault{10^{%d}}$' % i for i in range(-3, 1)],
229  [str(i) for i in range(1, 11)]]
230  xticks[0][-1] = ""
231  xticks[1][-1] = ""
232  for j in range(0, 4):
233  axes[1, j].set_xticklabels(xticks[j%2])
234  for t in axes[1, j].get_xticklabels():
235  t.set_fontsize(11)
236  fig.legend(handles, ["exp/dev", "lux/luv", "approx exp/dev", "approx lux/luv"],
237  loc='lower center', ncol=4)
238  fig.text(centers[0], 0.95, "exponential", ha='center', weight='bold')
239  fig.text(centers[1], 0.95, "de Vaucouleur", ha='center', weight='bold')
240  return fig, axes
241 
def evaluateRadial(basis, r, sbNormalize=False, doComponents=False)
Definition: tractor.py:92
def plotSuite(doComponents=False)
Definition: tractor.py:155

◆ registerRadialProfiles()

def lsst.shapelet.tractor.registerRadialProfiles ( )
Register the pickled profiles in the data directory with the RadialProfile
singleton registry.

This should only be called at import time by this module; it's only a function to
avoid polluting the module namespace with all the local variables used here.

Definition at line 47 of file tractor.py.

48  """Register the pickled profiles in the data directory with the RadialProfile
49  singleton registry.
50 
51  This should only be called at import time by this module; it's only a function to
52  avoid polluting the module namespace with all the local variables used here.
53  """
54  dataDir = os.path.join(os.environ["SHAPELET_DIR"], "data")
55  regex = re.compile(r"([a-z]+\d?)_K(\d+)_MR(\d+)\.pickle")
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)
66  warnings.warn("No C++ profile for multi-Gaussian pickle file '%s'" % filename)
67  continue
68  with open(os.path.join(dataDir, filename), 'rb') as stream:
69  if sys.version_info[0] >= 3:
70  array = pickle.load(stream, encoding='latin1')
71  else:
72  array = pickle.load(stream)
73  amplitudes = array[:nComponents]
74  amplitudes /= amplitudes.sum()
75  variances = array[nComponents:]
76  if amplitudes.shape != (nComponents,) or variances.shape != (nComponents,):
77  warnings.warn("Unknown format for multi-Gaussian pickle file '%s'" % filename)
78  continue
79  basis = MultiShapeletBasis(1)
80  for amplitude, variance in zip(amplitudes, variances):
81  radius = variance**0.5
82  matrix = numpy.array([[amplitude / ShapeletFunction.FLUX_FACTOR]], dtype=float)
83  basis.addComponent(radius, 0, matrix)
84  profile.registerBasis(basis, nComponents, maxRadius)
85 
86 
87 # We register all the profiles at module import time, to allow C++ code to access all available profiles
88 # without having to later call Python code to unpickle them.
90 
91 
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
def registerRadialProfiles()
Definition: tractor.py:47