LSSTApplications  18.0.0+106,18.0.0+50,19.0.0,19.0.0+1,19.0.0+10,19.0.0+11,19.0.0+13,19.0.0+17,19.0.0+2,19.0.0-1-g20d9b18+6,19.0.0-1-g425ff20,19.0.0-1-g5549ca4,19.0.0-1-g580fafe+6,19.0.0-1-g6fe20d0+1,19.0.0-1-g7011481+9,19.0.0-1-g8c57eb9+6,19.0.0-1-gb5175dc+11,19.0.0-1-gdc0e4a7+9,19.0.0-1-ge272bc4+6,19.0.0-1-ge3aa853,19.0.0-10-g448f008b,19.0.0-12-g6990b2c,19.0.0-2-g0d9f9cd+11,19.0.0-2-g3d9e4fb2+11,19.0.0-2-g5037de4,19.0.0-2-gb96a1c4+3,19.0.0-2-gd955cfd+15,19.0.0-3-g2d13df8,19.0.0-3-g6f3c7dc,19.0.0-4-g725f80e+11,19.0.0-4-ga671dab3b+1,19.0.0-4-gad373c5+3,19.0.0-5-ga2acb9c+2,19.0.0-5-gfe96e6c+2,w.2020.01
LSSTDataManagementBasePackage
densityPlot.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 """A set of matplotlib-based classes that displays a grid of 1-d and 2-d slices through an
24 N-d density.
25 
26 The main class, DensityPlot, manages the grid of matplotlib.axes.Axes objects, and holds
27 a sequence of Layer objects that each know how to draw individual 1-d or 2-d plots and a
28 data object that abstracts away how the N-d density data is actually represented.
29 
30 For simple cases, users can just create a custom data class with an interface like that of
31 the ExampleData class provided here, and use the provided HistogramLayer and SurfaceLayer
32 classes directly. In more complicated cases, users may want to create their own Layer classes,
33 which may define their own relationship with the data object.
34 """
35 
36 import collections.abc
37 import numpy
38 import matplotlib.cm
39 import matplotlib.pyplot
40 import matplotlib.ticker
41 
42 __all__ = ("HistogramLayer", "SurfaceLayer", "ScatterLayer", "CrossPointsLayer",
43  "DensityPlot", "ExampleData", "demo")
44 
45 
46 def hide_xticklabels(axes):
47  for label in axes.get_xticklabels():
48  label.set_visible(False)
49 
50 
51 def hide_yticklabels(axes):
52  for label in axes.get_yticklabels():
53  label.set_visible(False)
54 
55 
56 def mergeDefaults(kwds, defaults):
57  copy = defaults.copy()
58  if kwds is not None:
59  copy.update(**kwds)
60  return copy
61 
62 
64  """A Layer class for DensityPlot for gridded histograms, drawing bar plots in 1-d and
65  colormapped large-pixel images in 2-d.
66 
67  Relies on two data object attributes:
68 
69  values ----- a (M,N) array of data points, where N is the dimension of the dataset and M is the
70  number of data points
71 
72  weights ---- (optional) an array of weights with shape (M,); if not present, all weights will
73  be set to unity
74 
75  The need for these data object attributes can be removed by subclassing HistogramLayer and overriding
76  the hist1d and hist2d methods.
77  """
78 
79  defaults1d = dict(facecolor='b', alpha=0.5)
80  defaults2d = dict(cmap=matplotlib.cm.Blues, vmin=0.0, interpolation='nearest')
81 
82  def __init__(self, tag, bins1d=20, bins2d=(20, 20), kwds1d=None, kwds2d=None):
83  self.tag = tag
84  self.bins1d = bins1d
85  self.bins2d = bins2d
86  self.kwds1d = mergeDefaults(kwds1d, self.defaults1d)
87  self.kwds2d = mergeDefaults(kwds2d, self.defaults2d)
88 
89  def hist1d(self, data, dim, limits):
90  """Extract points from the data object and compute a 1-d histogram.
91 
92  Return value should match that of numpy.histogram: a tuple of (hist, edges),
93  where hist is a 1-d array with size=bins1d, and edges is a 1-d array with
94  size=self.bins1d+1 giving the upper and lower edges of the bins.
95  """
96  i = data.dimensions.index(dim)
97  if hasattr(data, "weights") and data.weights is not None:
98  weights = data.weights
99  else:
100  weights = None
101  return numpy.histogram(data.values[:, i], bins=self.bins1d, weights=weights,
102  range=limits, normed=True)
103 
104  def hist2d(self, data, xDim, yDim, xLimits, yLimits):
105  """Extract points from the data object and compute a 1-d histogram.
106 
107  Return value should match that of numpy.histogram2d: a tuple of (hist, xEdges, yEdges),
108  where hist is a 2-d array with shape=bins2d, xEdges is a 1-d array with size=bins2d[0]+1,
109  and yEdges is a 1-d array with size=bins2d[1]+1.
110  """
111  i = data.dimensions.index(yDim)
112  j = data.dimensions.index(xDim)
113  if hasattr(data, "weights") and data.weights is not None:
114  weights = data.weights
115  else:
116  weights = None
117  return numpy.histogram2d(data.values[:, j], data.values[:, i], bins=self.bins2d, weights=weights,
118  range=(xLimits, yLimits), normed=True)
119 
120  def plotX(self, axes, data, dim):
121  y, xEdge = self.hist1d(data, dim, axes.get_xlim())
122  xCenter = 0.5*(xEdge[:-1] + xEdge[1:])
123  width = xEdge[1:] - xEdge[:-1]
124  return axes.bar(xCenter, y, width=width, align='center', **self.kwds1d)
125 
126  def plotY(self, axes, data, dim):
127  x, yEdge = self.hist1d(data, dim, axes.get_ylim())
128  yCenter = 0.5*(yEdge[:-1] + yEdge[1:])
129  height = yEdge[1:] - yEdge[:-1]
130  return axes.barh(yCenter, x, height=height, align='center', **self.kwds1d)
131 
132  def plotXY(self, axes, data, xDim, yDim):
133  z, xEdge, yEdge = self.hist2d(data, xDim, yDim, axes.get_xlim(), axes.get_ylim())
134  return axes.imshow(z.transpose(), aspect='auto', extent=(xEdge[0], xEdge[-1], yEdge[0], yEdge[-1]),
135  origin='lower', **self.kwds2d)
136 
137 
139  """A Layer class that plots individual points in 2-d, and does nothing in 1-d.
140 
141  Relies on two data object attributes:
142 
143  values ----- a (M,N) array of data points, where N is the dimension of the dataset and M is the
144  number of data points
145 
146  weights ---- (optional) an array of weights with shape (M,); will be used to set the color of points
147 
148  """
149 
150  defaults = dict(linewidth=0, alpha=0.2)
151 
152  def __init__(self, tag, **kwds):
153  self.tag = tag
154  self.kwds = mergeDefaults(kwds, self.defaults)
155 
156  def plotX(self, axes, data, dim):
157  pass
158 
159  def plotY(self, axes, data, dim):
160  pass
161 
162  def plotXY(self, axes, data, xDim, yDim):
163  i = data.dimensions.index(yDim)
164  j = data.dimensions.index(xDim)
165  if hasattr(data, "weights") and data.weights is not None:
166  args = data.values[:, j], data.values[:, i], data.weights
167  else:
168  args = data.values[:, j], data.values[:, i]
169  return axes.scatter(*args, **self.kwds)
170 
171 
173  """A Layer class for analytic N-d distributions that can be evaluated in 1-d or 2-d slices.
174 
175  The 2-d slices are drawn as contours, and the 1-d slices are drawn as simple curves.
176 
177  Relies on eval1d and eval2d methods in the data object; this can be avoided by subclassing
178  SurfaceLayer and reimplementing its own eval1d and eval2d methods.
179  """
180 
181  defaults1d = dict(linewidth=2, color='r')
182  defaults2d = dict(linewidths=2, cmap=matplotlib.cm.Reds)
183 
184  def __init__(self, tag, steps1d=200, steps2d=200, filled=False, kwds1d=None, kwds2d=None):
185  self.tag = tag
186  self.steps1d = int(steps1d)
187  self.steps2d = int(steps2d)
188  self.filled = bool(filled)
189  self.kwds1d = mergeDefaults(kwds1d, self.defaults1d)
190  self.kwds2d = mergeDefaults(kwds2d, self.defaults2d)
191 
192  def eval1d(self, data, dim, x):
193  """Return analytic function values for the given values."""
194  return data.eval1d(dim, x)
195 
196  def eval2d(self, data, xDim, yDim, x, y):
197  """Return analytic function values for the given values."""
198  return data.eval2d(xDim, yDim, x, y)
199 
200  def plotX(self, axes, data, dim):
201  xMin, xMax = axes.get_xlim()
202  x = numpy.linspace(xMin, xMax, self.steps1d)
203  z = self.eval1d(data, dim, x)
204  if z is None:
205  return
206  return axes.plot(x, z, **self.kwds1d)
207 
208  def plotY(self, axes, data, dim):
209  yMin, yMax = axes.get_ylim()
210  y = numpy.linspace(yMin, yMax, self.steps1d)
211  z = self.eval1d(data, dim, y)
212  if z is None:
213  return
214  return axes.plot(z, y, **self.kwds1d)
215 
216  def plotXY(self, axes, data, xDim, yDim):
217  xMin, xMax = axes.get_xlim()
218  yMin, yMax = axes.get_ylim()
219  xc = numpy.linspace(xMin, xMax, self.steps2d)
220  yc = numpy.linspace(yMin, yMax, self.steps2d)
221  xg, yg = numpy.meshgrid(xc, yc)
222  z = self.eval2d(data, xDim, yDim, xg, yg)
223  if z is None:
224  return
225  if self.filled:
226  return axes.contourf(xg, yg, z, 6, **self.kwds2d)
227  else:
228  return axes.contour(xg, yg, z, 6, **self.kwds2d)
229 
230 
232  """A layer that marks a few points with axis-length vertical and horizontal lines.
233 
234  This relies on a "points" data object attribute.
235  """
236 
237  defaults = dict(alpha=0.8)
238 
239  def __init__(self, tag, colors=("y", "m", "c", "r", "g", "b"), **kwds):
240  self.tag = tag
241  self.colors = colors
242  self.kwds = mergeDefaults(kwds, self.defaults)
243 
244  def plotX(self, axes, data, dim):
245  i = data.dimensions.index(dim)
246  artists = []
247  for n, point in enumerate(data.points):
248  artists.append(axes.axvline(point[i], color=self.colors[n % len(self.colors)], **self.kwds))
249  return artists
250 
251  def plotY(self, axes, data, dim):
252  i = data.dimensions.index(dim)
253  artists = []
254  for n, point in enumerate(data.points):
255  artists.append(axes.axhline(point[i], color=self.colors[n % len(self.colors)], **self.kwds))
256  return artists
257 
258  def plotXY(self, axes, data, xDim, yDim):
259  i = data.dimensions.index(yDim)
260  j = data.dimensions.index(xDim)
261  artists = []
262  for n, point in enumerate(data.points):
263  artists.append(axes.axvline(point[j], color=self.colors[n % len(self.colors)], **self.kwds))
264  artists.append(axes.axhline(point[i], color=self.colors[n % len(self.colors)], **self.kwds))
265  return artists
266 
267 
269  """An object that manages a matrix of matplotlib.axes.Axes objects that represent a set of 1-d and 2-d
270  slices through an N-d density.
271  """
272 
273  class LayerDict(collections.abc.MutableMapping):
274 
275  def __init__(self, parent):
276  self._dict = dict()
277  self._parent = parent
278 
279  def __delitem__(self, name):
280  layer = self._dict.pop(name)
281  self._parent._dropLayer(name, layer)
282 
283  def __setitem__(self, name, layer):
284  self.pop(name, None)
285  self._dict[name] = layer
286  self._parent._plotLayer(name, layer)
287 
288  def __getitem__(self, name):
289  return self._dict[name]
290 
291  def __iter__(self):
292  return iter(self._dict)
293 
294  def __len__(self):
295  return len(self._dict)
296 
297  def __str__(self):
298  return str(self._dict)
299 
300  def __repr__(self):
301  return repr(self._dict)
302 
303  def replot(self, name):
304  layer = self._dict[name]
305  self._parent._dropLayer(name, layer)
306  self._parent._plotLayer(name, layer)
307 
308  def __init__(self, figure, **kwds):
309  self.figure = figure
310  self.data = dict(kwds)
311  active = []
312  self._lower = dict()
313  self._upper = dict()
314  # We merge the dimension name lists manually rather than using sets to preserve the order.
315  # Most of the time we expect all data objects to have the same dimensions anyway.
316  for v in self.data.values():
317  for dim in v.dimensions:
318  if dim not in active:
319  active.append(dim)
320  self._lower[dim] = v.lower[dim]
321  self._upper[dim] = v.upper[dim]
322  else:
323  self._lower[dim] = min(v.lower[dim], self._lower[dim])
324  self._upper[dim] = max(v.upper[dim], self._upper[dim])
325  self._active = tuple(active)
326  self._all_dims = frozenset(self._active)
327  self.figure.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95, hspace=0.01, wspace=0.01)
328  self._build_axes()
329  self.layers = self.LayerDict(self)
330 
331  def _dropLayer(self, name, layer):
332  def removeArtist(*key):
333  try:
334  self._objs.pop(key).remove()
335  except AttributeError:
336  # sometimes the value might be None, which doesn't have a remove
337  pass
338  except TypeError:
339  # probably a matplotlib bug: remove sometimes raises an exception,
340  # but it still works
341  pass
342  for i, yDim in enumerate(self._active):
343  removeArtist(None, i, name)
344  removeArtist(i, None, name)
345  for j, xDim in enumerate(self._active):
346  if i == j:
347  continue
348  removeArtist(i, j, name)
349 
350  def _plotLayer(self, name, layer):
351  for i, yDim in enumerate(self._active):
352  if yDim not in self.data[layer.tag].dimensions:
353  continue
354  self._objs[None, i, name] = layer.plotX(self._axes[None, i], self.data[layer.tag], yDim)
355  self._objs[i, None, name] = layer.plotY(self._axes[i, None], self.data[layer.tag], yDim)
356  for j, xDim in enumerate(self._active):
357  if xDim not in self.data[layer.tag].dimensions:
358  continue
359  if i == j:
360  continue
361  self._objs[i, j, name] = layer.plotXY(self._axes[i, j], self.data[layer.tag], xDim, yDim)
362  self._axes[None, i].xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(nbins=5, prune='both'))
363  self._axes[i, None].yaxis.set_major_locator(matplotlib.ticker.MaxNLocator(nbins=5, prune='both'))
364  self._axes[None, i].xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
365  self._axes[i, None].yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
366 
367  def _get_active(self):
368  return self._active
369 
370  def _set_active(self, active):
371  s = set(active)
372  if len(s) != len(active):
373  raise ValueError("Active set contains duplicates")
374  if not self._all_dims.issuperset(s):
375  raise ValueError("Invalid values in active set")
376  self._active = tuple(active)
377  self.replot()
378  active = property(_get_active, _set_active, doc="sequence of active dimensions to plot (sequence of str)")
379 
380  def replot(self):
381  self._lower = {dim: min(self.data[k].lower[dim] for k in self.data) for dim in self._active}
382  self._upper = {dim: max(self.data[k].upper[dim] for k in self.data) for dim in self._active}
383  self._build_axes()
384  for name, layer in self.layers.items():
385  self._plotLayer(name, layer)
386 
387  def _build_axes(self):
388  self.figure.clear()
389  self._axes = dict()
390  self._objs = dict()
391  n = len(self._active)
392  iStride = n + 1
393  jStride = -1
394  iStart = n + 1
395  jStart = n
396  for i in range(n):
397  j = i
398  axesX = self._axes[None, j] = self.figure.add_subplot(n+1, n+1, jStart+j*jStride)
399  axesX.autoscale(False, axis='x')
400  axesX.xaxis.tick_top()
401  axesX.set_xlim(self._lower[self._active[j]], self._upper[self._active[j]])
402  hide_yticklabels(axesX)
403  bbox = axesX.get_position()
404  bbox.y1 -= 0.035
405  axesX.set_position(bbox)
406  axesY = self._axes[i, None] = self.figure.add_subplot(n+1, n+1, iStart + iStart+i*iStride)
407  axesY.autoscale(False, axis='y')
408  axesY.yaxis.tick_right()
409  axesY.set_ylim(self._lower[self._active[i]], self._upper[self._active[i]])
410  hide_xticklabels(axesY)
411  bbox = axesY.get_position()
412  bbox.x1 -= 0.035
413  axesY.set_position(bbox)
414  for i in range(n):
415  for j in range(n):
416  axesXY = self._axes[i, j] = self.figure.add_subplot(
417  n+1, n+1, iStart+i*iStride + jStart+j*jStride,
418  sharex=self._axes[None, j],
419  sharey=self._axes[i, None]
420  )
421  axesXY.autoscale(False)
422  if j < n - 1:
423  hide_yticklabels(axesXY)
424  if i < n - 1:
425  hide_xticklabels(axesXY)
426  for i in range(n):
427  j = i
428  xbox = self._axes[None, j].get_position()
429  ybox = self._axes[i, None].get_position()
430  self.figure.text(0.5*(xbox.x0 + xbox.x1), 0.5*(ybox.y0 + ybox.y1), self.active[i],
431  ha='center', va='center', weight='bold')
432  self._axes[i, j].get_frame().set_facecolor('none')
433 
434  def draw(self):
435  self.figure.canvas.draw()
436 
437 
439  """An example data object for DensityPlot, demonstrating the necessarity interface.
440 
441  There are two levels of requirements for a data object. First are the attributes
442  required by the DensityPlot object itself; these must be present on every data object:
443 
444  dimensions ------ a sequence of strings that provide names for the dimensions
445 
446  lower ----------- a dictionary of {dimension-name: lower-bound}
447 
448  upper ----------- a dictionary of {dimension-name: upper-bound}
449 
450  The second level of requirements are those of the Layer objects provided here. These
451  may be absent if the associated Layer is not used or is subclassed to reimplement the
452  Layer method that calls the data object method. Currently, these include:
453 
454  eval1d, eval2d -- methods used by the SurfaceLayer class; see their docs for more info
455 
456  values ---------- attribute used by the HistogramLayer and ScatterLayer classes, an array
457  with shape (M,N), where N is the number of dimension and M is the number
458  of data points
459 
460  weights --------- optional attribute used by the HistogramLayer and ScatterLayer classes,
461  a 1-d array with size=M that provides weights for each data point
462  """
463 
464  def __init__(self):
465  self.dimensions = ["a", "b", "c"]
466  self.mu = numpy.array([-10.0, 0.0, 10.0])
467  self.sigma = numpy.array([3.0, 2.0, 1.0])
468  self.lower = {dim: -3*self.sigma[i] + self.mu[i] for i, dim in enumerate(self.dimensions)}
469  self.upper = {dim: 3*self.sigma[i] + self.mu[i] for i, dim in enumerate(self.dimensions)}
470  self.values = numpy.random.randn(2000, 3) * self.sigma[numpy.newaxis, :] + self.mu[numpy.newaxis, :]
471 
472  def eval1d(self, dim, x):
473  """Evaluate the 1-d analytic function for the given dim at points x (a 1-d numpy array;
474  this method must be numpy-vectorized).
475  """
476  i = self.dimensions.index(dim)
477  return numpy.exp(-0.5*((x-self.mu[i])/self.sigma[i])**2) / ((2.0*numpy.pi)**0.5 * self.sigma[i])
478 
479  def eval2d(self, xDim, yDim, x, y):
480  """Evaluate the 2-d analytic function for the given xDim and yDim at points x,y
481  (2-d numpy arrays with the same shape; this method must be numpy-vectorized).
482  """
483  i = self.dimensions.index(yDim)
484  j = self.dimensions.index(xDim)
485  return (numpy.exp(-0.5*(((x-self.mu[j])/self.sigma[j])**2 + ((y-self.mu[i])/self.sigma[i])**2)) /
486  (2.0*numpy.pi * self.sigma[j]*self.sigma[i]))
487 
488 
489 def demo():
490  """Create and return a DensityPlot with example data."""
491  fig = matplotlib.pyplot.figure()
492  p = DensityPlot(fig, primary=ExampleData())
493  p.layers['histogram'] = HistogramLayer('primary')
494  p.layers['surface'] = SurfaceLayer('primary')
495  p.draw()
496  return p
def hist2d(self, data, xDim, yDim, xLimits, yLimits)
Definition: densityPlot.py:104
def __init__(self, tag, colors=("y", "m", "c", "r", "g", "b"), kwds)
Definition: densityPlot.py:239
std::vector< SchemaItem< Flag > > * items
daf::base::PropertySet * set
Definition: fits.cc:902
def __init__(self, tag, bins1d=20, bins2d=(20, 20), kwds1d=None, kwds2d=None)
Definition: densityPlot.py:82
int min
int max
def eval2d(self, data, xDim, yDim, x, y)
Definition: densityPlot.py:196
def __init__(self, tag, steps1d=200, steps2d=200, filled=False, kwds1d=None, kwds2d=None)
Definition: densityPlot.py:184
def plotXY(self, axes, data, xDim, yDim)
Definition: densityPlot.py:216
def plotXY(self, axes, data, xDim, yDim)
Definition: densityPlot.py:162