#96 Binned WhiskerPlots

CHANGELOG.md

@@ -1,3 +1,4 @@
+7/13/18: Add BinnedScatterPlot systematics test type (issue #96)
 7/11/18: Update code to python3
 3/17/16: Update documentation to Sphinx standard and add documentation build files (issue #17)
 2/17/16: Changes to correlation function plots & documentation (issue #77)

stile/binning.py

@@ -23,10 +23,10 @@ class BinList:
     :returns:         A list of :class:`SingleBin` objects determined by the input criteria.
     """
 
-    def __init__(self, field, bin_list):
-        if not isinstance(field, str):
-            raise TypeError('Field description must be a string. Passed value: '+str(field)+
-                              'of type'+type(field))
+    def __init__(self, bin_list, field=None):
+        if not isinstance(field, str) and field is not None:
+            raise TypeError('Field description must be a string or None. Passed value: '+str(field)+
+                              ' of type'+str(type(field)))
         if not bin_list:
             raise TypeError('Must pass a non-empty bin_list')
         self.field = field
@@ -46,8 +46,12 @@ def __call__(self):
         Returns a list of :class:`SingleBin` objects following the definitions provided when the
         class was created.
         """
-        return_list = [SingleBin(field=self.field, low=low, high=high, short_name=str(i))
-                       for i, (low, high) in enumerate(zip(self.bin_list[:-1], self.bin_list[1:]))]
+        if self.field is not None:
+            bin_type = SingleBin
+        else:
+            bin_type = SingleBinFieldless
+        return_list = [bin_type(field=self.field, low=low, high=high, short_name=str(i))
+                       for i, (low, high) in  enumerate(zip(self.bin_list[:-1], self.bin_list[1:]))]
         if self.reverse:
             return_list.reverse()
         return return_list
@@ -79,10 +83,10 @@ class BinStep:
     :returns:         A list of :class:`SingleBin` objects determined by the input criteria.
     """
 
-    def __init__(self, field, low=None, high=None, step=None, n_bins=None, use_log=False):
+    def __init__(self, field=None, low=None, high=None, step=None, n_bins=None, use_log=False):
         if not isinstance(field, str):
             raise TypeError('Field description must be a string. Passed value: '+str(field)+
-                              'of type'+type(field))
+                              ' of type'+str(type(field)))
         self.field = field
         n_none = (low is None) + (high is None) + (step is None) + (n_bins is None)
         if n_none > 1:
@@ -143,20 +147,24 @@ def __init__(self, field, low=None, high=None, step=None, n_bins=None, use_log=F
             self.reverse = False
 
     def __call__(self):
+        if self.field is not None:
+            bin_type = SingleBin
+        else:
+            bin_type = SingleBinFieldless
         if self.use_log:
-            return_list = [SingleBin(field=self.field, low=numpy.exp(self.low+i*self.step),
+            return_list = [bin_type(field=self.field, low=numpy.exp(self.low+i*self.step),
                                      high=numpy.exp(self.low+(i+1)*self.step),
                                      short_name=str(i)) for i in range(self.n_bins)]
         else:
-            return_list = [SingleBin(field=self.field, low=self.low+i*self.step,
+            return_list = [bin_type(field=self.field, low=self.low+i*self.step,
                                      high=self.low+(i+1)*self.step, short_name=str(i))
                                      for i in range(self.n_bins)]
         if self.reverse:
             return_list.reverse()
         return return_list
 
 
-class SingleBin:
+class SingleBinFieldless(object):
     """
     A class that contains the information for one particular bin generated from one of
     the :class:`Bin*` classes. The attributes can be accessed directly for :class:`DataHandler`\s
@@ -173,15 +181,11 @@ class SingleBin:
                        [default: ``"low-high"``].
     """
 
-    def __init__(self, field, low, high, short_name, long_name=None):
-        if not isinstance(field, str):
-            raise TypeError('Field description must be a string. Passed value: '+str(field)+
-                              'of type'+type(field))
+    def __init__(self, low, high, short_name, long_name=None):
         if high <= low:
             raise ValueError("High ("+str(high)+") must be greater than low ("+str(low)+")")
         if not isinstance(short_name, str) or (long_name and not isinstance(long_name, str)):
             raise TypeError("Short_name and long_name must be strings")
-        self.field = field
         self.low = low
         self.high = high
         self.short_name = short_name
@@ -199,7 +203,46 @@ def __call__(self, data):
         :returns:      A NumPy array corresponding to the input data, restricted to the bin
                        described by this object.
         """
-        return data[numpy.logical_and(data[self.field] >= self.low, data[self.field] < self.high)]
+        return data[numpy.logical_and(data >= self.low, data < self.high)]
+
+
+class SingleBin(SingleBinFieldless):
+    """
+    A class that contains the information for one particular bin generated from one of the Bin*
+    classes. The attributes can be accessed directly for DataHandlers that read in the data
+    selectively. The class can also be called with a data array to bin it to the correct data
+    range: SingleBin(array) will return only the data within the bounds of the particular instance
+    of the class.  The endpoints are assumed to be [low,high), that is, low <= data < high, with
+    defined relational operators.
+
+    This version is for data arrays with defined fields such as numpy.recarrays; see 
+    :class:`SingleBinFieldless` for the version that operates on raw arrays.
+
+    @param field      The index of the field containing the data to be binned (must be a string).
+    @param low        The lower edge of the bin (inclusive).
+    @param high       The upper edge of the bin (exclusive).
+    @param short_name A string denoting this bin in filenames.
+    @param long_name  A string denoting this bin in program text outputs/plots.
+                      [default: "low-high"]
+    """
+    def __init__(self, field, low, high, short_name, long_name=None):
+        if not isinstance(field, str):
+            raise TypeError('Field description must be a string. Passed value: '+str(field)+
+                              ' of type'+str(type(field)))
+        self.field = field
+        super(SingleBin, self).__init__(low, high, short_name, long_name)
+
+    def __call__(self, data):
+        """
+        Given data, returns only the data with data[self.field] within the bounds
+        [self.low, self.high).
+
+        @param data   A NumPy array of data that can be indexed by self.field.
+        @returns      A NumPy array corresponding to the input data, restricted to the bin
+                      described by this object.
+        """
+        return super(SingleBin, self).__call__(data[self.field])
+
 
 
 class BinFunction:

stile/sys_tests.py

@@ -1191,6 +1191,10 @@ def WhiskerPlotSysTest(type=None):
         - **Star**: whisker plot of shapes of PSF stars
         - **PSF**: whisker plot of PSF shapes at the location of PSF stars
         - **Residual**: whisker plot of (star shape-PSF shape)
+        - **BinnedStar**: whisker plot of shapes of PSF stars, based on averages in pixels
+        - **BinnedPSF**: whisker plot of PSF shapes at the location of PSF stars, based on averages in 
+          pixels
+        - **BinnedResidual**: whisker plot of (star shape-PSF shape), based on averages in pixels
         - **None**: an empty :class:`BaseWhiskerPlotSysTest` class instance, which can be used for
           multiple types of whisker plots.  See the documentation
           for :class:`BaseWhiskerPlotSysTest` (especially the method
@@ -1204,6 +1208,12 @@ def WhiskerPlotSysTest(type=None):
         return WhiskerPlotPSFSysTest()
     elif type=='Residual':
         return WhiskerPlotResidualSysTest()
+    elif type=='BinnedStar':
+        return BinnedWhiskerPlotStarSysTest()
+    elif type=='BinnedPSF':
+        return BinnedWhiskerPlotPSFSysTest()
+    elif type=='BinnedResidual':
+        return BinnedWhiskerPlotResidualSysTest()
     elif type is None:
         return BaseWhiskerPlotSysTest()
     else:
@@ -1321,7 +1331,7 @@ class WhiskerPlotStarSysTest(BaseWhiskerPlotSysTest):
     required_quantities = [('x', 'y', 'g1', 'g2', 'sigma')]
 
     def __call__(self, array, linewidth=0.01, scale=None, figsize=None,
-                 xlim=None, ylim=None):
+                 xlim=None, ylim=None, **kwargs):
         if 'CCD' in array.dtype.names:
             fields = list(self.required_quantities[0]) + ['CCD']
         else:
@@ -1331,7 +1341,7 @@ def __call__(self, array, linewidth=0.01, scale=None, figsize=None,
                                 linewidth=linewidth, scale=scale, figsize=figsize,
                                 xlabel=r'$x$ [pixel]', ylabel=r'$y$ [pixel]',
                                 size_label=r'$\sigma$ [pixel]',
-                                xlim=xlim, ylim=ylim, equal_axis=True)
+                                xlim=xlim, ylim=ylim, equal_axis=True, **kwargs)
 
 
 class WhiskerPlotPSFSysTest(BaseWhiskerPlotSysTest):
@@ -1344,7 +1354,7 @@ class WhiskerPlotPSFSysTest(BaseWhiskerPlotSysTest):
     required_quantities = [('x', 'y', 'psf_g1', 'psf_g2', 'psf_sigma')]
 
     def __call__(self, array, linewidth=0.01, scale=None, figsize=None,
-                 xlim=None, ylim=None):
+                 xlim=None, ylim=None, **kwargs):
         if 'CCD' in array.dtype.names:
             fields = list(self.required_quantities[0]) + ['CCD']
         else:
@@ -1354,7 +1364,7 @@ def __call__(self, array, linewidth=0.01, scale=None, figsize=None,
                                 array['psf_sigma'], linewidth=linewidth, scale=scale,
                                 figsize=figsize, xlabel=r'$x$ [pixel]', ylabel=r'$y$ [pixel]',
                                 size_label=r'$\sigma$ [pixel]',
-                                xlim=xlim, ylim=ylim, equal_axis=True)
+                                xlim=xlim, ylim=ylim, equal_axis=True, **kwargs)
 
 
 class WhiskerPlotResidualSysTest(BaseWhiskerPlotSysTest):
@@ -1366,7 +1376,7 @@ class WhiskerPlotResidualSysTest(BaseWhiskerPlotSysTest):
     required_quantities = [('x', 'y', 'g1', 'g2', 'sigma', 'psf_g1', 'psf_g2', 'psf_sigma')]
 
     def __call__(self, array, linewidth=0.01, scale=None, figsize=None,
-                 xlim=None, ylim=None):
+                 xlim=None, ylim=None, **kwargs):
         data = [array['x'], array['y'], array['g1'] - array['psf_g1'],
                 array['g2'] - array['psf_g2'], array['sigma'] - array['psf_sigma']]
         fields = ['x', 'y', 'g1-psf_g1', 'g2-psf_g2', 'sigma-psf_sigma']
@@ -1379,7 +1389,71 @@ def __call__(self, array, linewidth=0.01, scale=None, figsize=None,
                                 linewidth=linewidth, scale=scale,
                                 figsize=figsize, xlabel=r'$x$ [pixel]', ylabel=r'$y$ [pixel]',
                                 size_label=r'$\sigma$ [pixel]',
-                                xlim=xlim, ylim=ylim, equal_axis=True)
+                                xlim=xlim, ylim=ylim, equal_axis=True, **kwargs)
+
+class BinnedWhiskerPlotSysTest(BaseWhiskerPlotSysTest):
+    def whiskerPlot(self, x, y, g1, g2, sigma, n_x=30, n_y=30, split_by_objects=True, 
+                                **kwargs):
+        """
+        A method for making whisker plots where the whiskers represent the average value in pixels.
+        Most of the arguments and keywords are the same as for 
+        :func:`BaseWhiskerPlotSysTest.whiskerPlot`, with the following additional keywords:
+
+        @param n_x              Number of pixels to use in the x-coordinate [default: 30]
+        @param n_y              Number of pixels to use in the y-coordinate [default: 30]
+        @param split_by_objects If True, make the pixels so there are an equal number of objects
+                                in each column or row. (Note this doesn't guarantee equal number
+                                in each pixel, just each column or row of pixels in total.)  If 
+                                False, instead use equal binning in the x- or y-coordinate, eg
+                                (0, 10, 20, 30...) for the bin endpoints.  [default: True]
+        """
+        new_x = []
+        new_y = []
+        new_g1 = []
+        new_g2 = []
+        new_sigma = []
+        if split_by_objects:
+            bin_x = numpy.percentile(x, numpy.linspace(0, 100, num=n_x+1, endpoint=True))
+            bin_y = numpy.percentile(y, numpy.linspace(0, 100, num=n_y+1, endpoint=True))
+        else:
+            bin_x = numpy.linspace(numpy.min(x), numpy.max(x), num=n_x, endpoint=True)
+            bin_y = numpy.linspace(numpy.min(y), numpy.max(y), num=n_y, endpoint=True)
+        for lo_x, hi_x in zip(bin_x[:-1], bin_x[1:]):
+            for lo_y, hi_y in zip(bin_y[:-1], bin_y[1:]):
+                mask = (x>lo_x) & (x<=hi_x) & (y>lo_y) & (y<=hi_y)
+                if numpy.any(mask): # guard against objectless pixells
+                    new_x.append(numpy.mean(x[mask]))
+                    new_y.append(numpy.mean(y[mask]))
+                    new_g1.append(numpy.mean(g1[mask]))
+                    new_g2.append(numpy.mean(g2[mask]))
+                    new_sigma.append(numpy.mean(sigma[mask]))
+        new_x = numpy.asarray(new_x)
+        new_y = numpy.asarray(new_y)
+        new_g1 = numpy.asarray(new_g1)
+        new_g2 = numpy.asarray(new_g2)
+        new_sigma = numpy.asarray(new_sigma)
+        return super(BinnedWhiskerPlotSysTest, self).whiskerPlot(
+                                        new_x, new_y, new_g1, new_g2, new_sigma, **kwargs)   
+
+
+class BinnedWhiskerPlotStarSysTest(BinnedWhiskerPlotSysTest, WhiskerPlotStarSysTest):
+    short_name = 'binned_whiskerplot_star'
+    long_name = 'Make a binned whisker plot of stars'
+    objects_list = ['star PSF']
+    required_quantities = [('x', 'y', 'g1', 'g2', 'sigma')]
+
+class BinnedWhiskerPlotPSFSysTest(BinnedWhiskerPlotSysTest, WhiskerPlotPSFSysTest):
+    short_name = 'binned_whiskerplot_psf'
+    long_name = 'Make a binned whisker plot of PSFs'
+    objects_list = ['star PSF']
+    required_quantities = [('x', 'y', 'psf_g1', 'psf_g2', 'psf_sigma')]
+
+class BinnedWhiskerPlotResidualSysTest(BinnedWhiskerPlotSysTest, WhiskerPlotResidualSysTest):
+    short_name = 'whiskerplot_residual'
+    long_name = 'Make a Whisker plot of residuals'
+    objects_list = ['star PSF']
+    required_quantities = [('x', 'y', 'g1', 'g2', 'sigma', 'psf_g1', 'psf_g2', 'psf_sigma')]
+
 
 class HistogramSysTest(SysTest):
     """