LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
interactive.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 
23 import numpy
24 
25 import lsst.pipe.base
26 import lsst.pex.config
28 import lsst.afw.image
31 
32 from ..measureCcd import MeasureCcdTask
33 from ..measureCoadd import MeasureCoaddTask
34 from ..measureMulti import MeasureMultiTask
35 from .. import modelfitLib
36 
37 from .densityPlot import CrossPointsLayer, DensityPlot, HistogramLayer, SurfaceLayer
38 from .modelFitAdapters import OptimizerDataAdapter, OptimizerTrackLayer, SamplingDataAdapter
39 from .optimizerDisplay import OptimizerDisplay
40 
41 __all__ = ("Interactive", )
42 
43 
45  """Interactive analysis helper class
46 
47  This class manages a butler, calexp, modelfits catalog, and an instance
48  of a Measure*Task, allowing individual objects to be re-fit and plotted.
49  """
50 
51  def __init__(self, root, tag=None, config=None, dataId=None, mode="ccd"):
52  """Construct an interactive analysis object.
53 
54  @param[in] rerun Output directory, relative to $S13_DATA_DIR/output.
55  measureCcd.py (or measureCoadd.py if mode='coadd') must
56  have been run (possibly with prepOnly=True) previously
57  with this output directory.
58  @param[in] tag Tag associated with the run; see BaseMeasureConfig.tag.
59  If None, config must not be (and config.tag will be used).
60  @param[in] config ConfigClass instance; if None, it will be loaded from disk.
61  @param[in] dataId Butler data ID of the image to analyze.
62  @param[in] mode One of "ccd", "coadd", or "multi", indicating whether
63  to use MeasureCcdTask, MeasureCoaddTask, or MeasureMultiTask.
64  """
65  self.modemode = mode.lower()
67  TaskClass = None
68  configName = None
69  if self.modemode == "ccd":
70  if dataId is None:
71  dataId = dict(visit=100, raft="2,2", sensor="1,1")
72  self.dataRefdataRef = self.butlerbutler.dataRef("calexp", dataId=dataId)
73  configName = "measureCcd_config"
74  TaskClass = MeasureCcdTask
75  elif self.modemode == "coadd":
76  if dataId is None:
77  dataId = dict(tract=0, patch="2,2", filter="i")
78  self.dataRefdataRef = self.butlerbutler.dataRef("deepCoadd_calexp", dataId=dataId)
79  configName = "deep_measureCoadd_config"
80  TaskClass = MeasureCoaddTask
81  elif self.modemode.startswith("multi"):
82  if dataId is None:
83  dataId = dict(tract=0, patch="2,2", filter="i")
84  self.dataRefdataRef = self.butlerbutler.dataRef("deepCoadd_calexp", dataId=dataId)
85  configName = "deep_measureMulti_config"
86  TaskClass = MeasureMultiTask
87  if config is None:
88  config = self.butlerbutler.get(configName, tag=tag, immediate=True)
89  config.previous = tag
90  if tag is None:
91  if config is None:
92  raise ValueError("tag and config arguments cannot both be None")
93  tag = config.tag
94  else:
95  config.tag = "intermediate"
96  self.tasktask = TaskClass(config=config, butler=self.butlerbutler)
97  self.inputsinputs = self.tasktask.readInputs(self.dataRefdataRef)
98 
99  def fit(self, outRecord):
100  """Re-fit the object indicated by the given record sequential index or source ID,
101  returning the record.
102  """
103  likelihood = self.tasktask.makeLikelihood(self.inputsinputs, outRecord)
104  self.tasktask.fitter.run(likelihood, outRecord)
105  return outRecord
106 
107  def plotDistribution(self, *records):
108  """Plot a representation of the posterior distribution from a ModelFitRecord.
109  """
110  import matplotlib
111  recordId = records[0].getId()
112  figure = matplotlib.pyplot.figure(recordId, figsize=(10, 10))
113  data = {}
114  layers = {}
115  for record in records:
116  assert record.getId() == recordId
117  if modelfitLib.MarginalSamplingInterpreter.cast(record.getInterpreter()):
118  data["marginal"] = SamplingDataAdapter(record)
119  layers["marginal.samples"] = HistogramLayer("direct")
120  layers["marginal.proposal"] = SurfaceLayer("direct")
121  elif modelfitLib.DirectSamplingInterpreter.cast(record.getInterpreter()):
122  data["direct"] = SamplingDataAdapter(record)
123  layers["direct.samples"] = HistogramLayer("direct")
124  layers["direct.proposal"] = SurfaceLayer("direct")
125  elif modelfitLib.OptimizerInterpreter.cast(record.getInterpreter()):
126  data["optimizer"] = OptimizerDataAdapter(record)
127  layers["optimizer.track"] = OptimizerTrackLayer("optimizer")
128  layers["optimizer.pdf"] = SurfaceLayer("optimizer",
129  kwds1d=dict(color='g'),
130  kwds2d=dict(cmap=matplotlib.cm.Greens))
131  layers["optimizer.points"] = CrossPointsLayer("optimizer")
132  else:
133  raise ValueError("Unknown or missing interpreter")
134  p = DensityPlot(figure, **data)
135  p.layers.update(layers)
136  p.draw()
137  return p
138 
139  def displayOptimizer(self, record, **kwds):
140  likelihood = self.tasktask.makeLikelihood(self.inputsinputs, record)
141  objective = modelfitLib.OptimizerObjective.makeFromLikelihood(likelihood, self.tasktask.prior)
142  return OptimizerDisplay(record.getSamples(), self.tasktask.model, objective)
143 
144  def displayResiduals(self, record, nonlinear="fit", amplitudes="fit", doApplyWeights=False):
145  """Display the data postage stamp along with the model image and residuals in ds9.
146 
147  @param[in] record ModelFitRecord defining the object to display
148  @param[in] nonlinear Vector of nonlinear parameters, or a string prefix (see below)
149  @param[in] amplitudes Vector of linear parameters, or a string prefix (see below)
150  @param[in] doApplyWeights Whether to show the weighted images used directly in the fit
151  or the original unweighted data.
152 
153  String prefixes are used to extract the parameters from the record. Usually the following
154  are available:
155  fit ------- use record.get("fit.*"); the best-fit parameters
156  initial --- use record.get("initial.*"); the initial parameters
157  """
158  likelihood = self.tasktask.makeLikelihood(self.inputsinputs, record)
159 
160  if isinstance(nonlinear, str):
161  nonlinear = record.get(nonlinear + ".nonlinear")
162  else:
163  assert nonlinear.shape == (likelihood.getNonlinearDim(),)
164 
165  matrix = numpy.zeros((likelihood.getAmplitudeDim(), likelihood.getDataDim()),
166  dtype=modelfitLib.Pixel).transpose()
167  likelihood.computeModelMatrix(matrix, nonlinear, doApplyWeights)
168 
169  if isinstance(amplitudes, str):
170  amplitudes = record.get(amplitudes + ".amplitudes")
171  else:
172  assert amplitudes.shape == (likelihood.getAmplitudeDim(),)
173 
174  bbox = record.getFootprint().getBBox()
175  bbox.grow(2)
176  flatModel = numpy.zeros(likelihood.getDataDim(), dtype=modelfitLib.Pixel)
177  flatModel[:] = numpy.dot(matrix, amplitudes)
178 
179  imgData = lsst.afw.image.MaskedImageF(self.inputsinputs.exposure.getMaskedImage(), bbox,
180  lsst.afw.image.PARENT, True)
181  bitmask = imgData.getMask().addMaskPlane("FIT_REGION")
182  regionMask = lsst.afw.image.MaskU(bbox)
183  lsst.afw.detection.setMaskFromFootprint(regionMask, record.getFootprint(), bitmask)
184  dataMask = imgData.getMask()
185  dataMask |= regionMask
186  if doApplyWeights:
187  imgData.getImage().set(0.0)
188  imgData.getVariance().set(0.0)
189  flatData = likelihood.getData()
190  lsst.afw.detection.expandArray(record.getFootprint(), flatData, imgData.getImage().getArray(),
191  imgData.getXY0())
192  imgModel = lsst.afw.image.MaskedImageF(lsst.afw.image.ImageF(bbox), regionMask)
193  lsst.afw.detection.expandArray(record.getFootprint(), flatModel, imgModel.getImage().getArray(),
194  imgModel.getXY0())
195  imgResiduals = lsst.afw.image.MaskedImageF(imgData, True)
196  imgResiduals -= imgModel
198  mosaic.setMode("x")
199  mosaic.append(imgData, "data")
200  mosaic.append(imgModel, "model")
201  mosaic.append(imgResiduals, "data-model")
202  grid = mosaic.makeMosaic()
def displayResiduals(self, record, nonlinear="fit", amplitudes="fit", doApplyWeights=False)
Definition: interactive.py:144
def __init__(self, root, tag=None, config=None, dataId=None, mode="ccd")
Definition: interactive.py:51
daf::base::PropertySet * set
Definition: fits.cc:912
def setMaskTransparency(name, frame=None)
Definition: ds9.py:80
def mtv(data, frame=None, title="", wcs=None, *args, **kwargs)
Definition: ds9.py:92
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.