diff --git a/CHANGELOG.md b/CHANGELOG.md index d07c664..58c574f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,4 @@ +7/12/18: Add BinnedScatterPlotSysTest class (issue #68) 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) diff --git a/stile/__init__.py b/stile/__init__.py index 22ab232..80d09bb 100644 --- a/stile/__init__.py +++ b/stile/__init__.py @@ -6,4 +6,5 @@ from .treecorr_utils import ReadTreeCorrResultsFile from .data_handler import DataHandler from .sys_tests import (StatSysTest, CorrelationFunctionSysTest, ScatterPlotSysTest, - WhiskerPlotSysTest, HistogramSysTest) + WhiskerPlotSysTest, HistogramSysTest, BinnedScatterPlotSysTest) + diff --git a/stile/hsc/base_tasks.py b/stile/hsc/base_tasks.py index ee06258..1cf1954 100644 --- a/stile/hsc/base_tasks.py +++ b/stile/hsc/base_tasks.py @@ -75,7 +75,8 @@ class CCDSingleEpochStileConfig(lsst.pex.config.Config): "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", - "ScatterPlotResidualSigmaVsPSFMag" + "ScatterPlotResidualSigmaVsPSFMag", + "RMSESky", "RMSEChip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control TreeCorr", keytype=str, itemtype=str, @@ -744,7 +745,8 @@ class VisitSingleEpochStileConfig(CCDSingleEpochStileConfig): "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", - "ScatterPlotResidualSigmaVsPSFMag" + "ScatterPlotResidualSigmaVsPSFMag", + "RMSESky", "RMSEChip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control treecorr", keytype=str, itemtype=str, @@ -1062,7 +1064,8 @@ class PatchSingleEpochStileConfig(CCDSingleEpochStileConfig): "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", - "ScatterPlotResidualSigmaVsPSFMag" + "ScatterPlotResidualSigmaVsPSFMag", + "RMSESky", "RMSEChip", "CountPerMagnitude" ]) do_hsm = lsst.pex.config.Field(dtype=bool, default=True, doc="Use HSM shapes for galaxies?") treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control TreeCorr", @@ -1159,7 +1162,8 @@ class TractSingleEpochStileConfig(CCDSingleEpochStileConfig): "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", - "ScatterPlotResidualSigmaVsPSFMag" + "ScatterPlotResidualSigmaVsPSFMag", + "RMSESky", "RMSEChip", "CountPerMagnitude" ]) do_hsm = lsst.pex.config.Field(dtype=bool, default=True, doc="Use HSM shapes for galaxies?") treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control treecorr", diff --git a/stile/hsc/sys_test_adapters.py b/stile/hsc/sys_test_adapters.py index 796dfc6..bc6710c 100644 --- a/stile/hsc/sys_test_adapters.py +++ b/stile/hsc/sys_test_adapters.py @@ -1,4 +1,4 @@ -""" +""" sys_test_adapters.py: Contains classes to wrap Stile systematics tests with functions necessary to run the tests via the HSC/LSST pipeline. """ @@ -474,6 +474,79 @@ def __call__(self, task_config, *data, **kwargs): new_data = [self.fixArray(d) for d in data] return self.sys_test(*new_data, per_ccd_stat=per_ccd_stat) +class RMSE1SkyAdapter(NoConfigShapeSysTestAdapter): + def __init__(self, config): + self.shape_type = 'sky' + self.config = config + self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', y_field='g1', w_field='w', method='rms') + self.name = 'rms_e1_sky' + self.setupMasks(['galaxy']) + +class RMSE2SkyAdapter(NoConfigShapeSysTestAdapter): + def __init__(self, config): + self.shape_type = 'sky' + self.config = config + self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', y_field='g2', w_field='w', method='rms') + self.name = 'rms_e2_sky' + self.setupMasks(['galaxy']) + +class RMSE1ChipAdapter(NoConfigShapeSysTestAdapter): + def __init__(self, config): + self.shape_type = 'chip' + self.config = config + self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', y_field='g1', w_field='w', method='rms') + self.name = 'rms_e1_chip' + self.setupMasks(['galaxy']) + +class RMSE2ChipAdapter(NoConfigShapeSysTestAdapter): + def __init__(self, config): + self.shape_type = 'chip' + self.config = config + self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', y_field='g2', w_field='w', method='rms') + self.name = 'rms_e2_chip' + self.setupMasks(['galaxy']) + +class CountPerMagnitudeAdapter(BaseSysTestAdapter): + def __init__(self, config): + self.config = config + self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', w_field='w', method='count') + self.name = 'count_per_mag' + self.setupMasks(['galaxy']) + +class RMSESkyAdapter(ShapeSysTestAdapter): + def __init__(self, config): + self.shape_type = 'sky' + self.config = config + self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', y_field='g', w_field='w', method='rms') + self.name = 'rms_e_sky' + self.setupMasks(['galaxy']) + + def __call__(self, task_config, *data, **kwargs): + new_data = [] + for d in data: + if d is None: + new_data.append(d) + else: + # We'll just make a new array using numpy.rec.fromarrays that includes |g|. + # First collect everything else, with a placeholder [] for |g| + list_of_arrays = [d[field] if field is not 'g' else [] for field in self.sys_test.required_quantities] + shape_col_1 = 'g1_'+self.shape_type + shape_col_2 = 'g2_'+self.shape_type + # Then replace that empty list with |g| + list_of_arrays[self.sys_test.required_quantities.index('g')] = ( + numpy.sqrt(d[shape_col_1]**2+d[shape_col_2]**2)) + new_data.append(numpy.rec.fromarrays(list_of_arrays, self.sys_test.required_quantities)) + return self.sys_test(*new_data, **kwargs) + +class RMSEChipAdapter(RMSESkyAdapter): + # We can use the same __call__ as RMSESKY, we just need to switch self.shape_type to 'chip' + def __init__(self, config): + self.shape_type = 'chip' + self.config = config + self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', y_field='g', w_field='w', method='rms') + self.name = 'rms_e_sky' + self.setupMasks(['galaxy']) + adapter_registry.register("StatsPSFFlux", StatsPSFFluxAdapter) adapter_registry.register("GalaxyShear", GalaxyShearAdapter) adapter_registry.register("BrightStarShear", BrightStarShearAdapter) @@ -492,3 +565,10 @@ def __call__(self, task_config, *data, **kwargs): adapter_registry.register("ScatterPlotResidualVsPSFSigma", ScatterPlotResidualVsPSFSigmaAdapter) adapter_registry.register("ScatterPlotResidualSigmaVsPSFMag", ScatterPlotResidualSigmaVsPSFMagAdapter) +adapter_registry.register("RMSE1Sky", RMSE1SkyAdapter) +adapter_registry.register("RMSE2Sky", RMSE2SkyAdapter) +adapter_registry.register("RMSE1Chip", RMSE1ChipAdapter) +adapter_registry.register("RMSE2Chip", RMSE2ChipAdapter) +adapter_registry.register("CountPerMagnitude", CountPerMagnitudeAdapter) +adapter_registry.register("RMSESky", RMSE1SkyAdapter) +adapter_registry.register("RMSEChip", RMSE1ChipAdapter) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index 276dbb9..1197323 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -2369,3 +2369,266 @@ def __call__(self, array, per_ccd_stat='None', color='', lim=None): color=color, lim=lim, equal_axis=False, linear_regression=True, reference_line='zero') +class BinnedScatterPlotSysTest(BaseScatterPlotSysTest): + short_name = 'binnedscatterplot' + """ + A base class for Stile systematics tests that generate scatter plots binned in one variable, + such as RMS ellipticity in magnitude bins or number density in longitude bins. This implements + a `__call__` method to produce a `matplotlib.figure.Figure`. Every child class of + BinnedScatterPlotSysTest should use binnedScatterPlotSysTest.scatterPlot through `super`. See + the docstring for BinnedScatterPlotSysTest.binnedScatterPlot for information on how to write + further tests using it. + + We automatically implement four statistics to report per bin: 'mean', 'median', 'rms' and + 'count'. 'mean' reports the mean of all y values in each bin; 'median' reports the median, with + an error estimate that is valid only for large numbers of points; 'rms' reports the root-mean- + square value, with an approximate error estimate; and 'count' simply report the number of + objects, with sqrt(n) errors. The user may also pass a callable function that operates on the + entire data array (not just the fields defined by `x_field` and `y_field`, described below, + which are used by the four methods described here). The function should return either a single + value that will be used for the bin, or a tuple of two numbers, the value for the bin and its + error. + + Weights or y errors can be given for any of the builtin methods except 'median'. Y errors are + assumed to be Gaussian and the derived weights are 1/yerr**2. Weighted `count`s are simply + the sums of the weights. + + The following quantities can be defined when the class is initialized, or passed to the + `__call__` method as kwargs: + @param x_field The name of the field in `array` to be used for x (the field that + defines the bins--magnitude in a plot of RMS ellipticity as a + function of magnitude, for example). + @param y_field The name of the field in `array` to be used for y. + @param yerr_field The name of the field in `array` to be used for y error. + @param w_field The name of the field in `array` to be used for the weights. + Defining both yerr_field and w_field will result in an error. + @param method The method to combine the values of `array[y_field]` in each bin. + Can be a string ('mean', 'median', 'rms' or 'count'), or a callable + function which operates on the given array and returns either the + desired value for the (already-binned) data, or a tuple of the value + and its error. [default: 'mean'] + @param binning The way to bin the values based on `array[x_field]`. If not given, + ten equally-sized bins will be used; if given, should be a number + indicating the number of requested equally-sized bins (will be + rounded to nearest int) or one of the Binning* classes defined in + Stile (BinStep, BinList, etc). + [default: None, meaning ten equally-sized bins] + """ + def __init__(self, x_field=None, y_field=None, yerr_field=None, w_field=None, method='mean', + binning=10): + self.x_field = x_field + self.y_field = y_field + self.yerr_field = yerr_field + self.w_field = w_field + list_of_quantities = [] + for quantity in [x_field, y_field, yerr_field, w_field]: + if quantity: + list_of_quantities.append(quantity) + self.required_quantities = [list_of_quantities] # need a list per object type + self.method = method + self.binning = binning + + def __call__(self, array, x_field=None, y_field=None, yerr_field=None, w_field=None, + method=None, binning=None, xlabel=None, ylabel=None, zlabel=None, color = "", + lim=None, equal_axis=False, linear_regression=False, reference_line = None, + area=None, area_units=None): + """ + Draw a binned scatter plot and return a `matplotlib.figure.Figure` object. + This method has a bunch of options for controlling appearance of a plot, which are + explained below. To implement a child class of BinnedScatterPlotSysTest, call + `super(classname, self).__call__(*args, **kwargs)` within the `__call__` method of the child + class and return the `matplotlib.figure.Figure` that this `__call__` method returns. + + Kwargs except `array`, `*_field`, `method` and `binning` are passed through to the parent + class function `scatterPlot()` unaltered. + + @param array A structured NumPy array which contains data to be plotted. + @param x_field The name of the field in `array` to be used for x (the field that + defines the bins--magnitude in a plot of RMS ellipticity as a + function of magnitude, for example). + @param y_field The name of the field in `array` to be used for y. + @param yerr_field The name of the field in `array` to be used for y error. + @param w_field The name of the field in `array` to be used for the weights. + Defining both yerr_field and w_field will result in an error. + @param method The method to combine the values of `array[y_field]` in each bin. + Can be a string ('mean', 'median', 'rms' or 'count'), or a callable + function which operates on the given array and returns either the + desired value for the (already-binned) data, or a tuple of the value + and its error. If 'count' is passed, the user can also pass the + kwarg 'area' to instead return a number density, and optionally also + a kwarg 'area_unit' for use in axis and column labels. + [default: 'mean'] + @param binning The way to bin the values based on `array[x_field]`. If not given, + ten equally-sized bins will be used; if given, should be a number + indicating the number of requested equally-sized bins (will be + rounded to nearest int) or one of the Binning* classes defined in + Stile (BinStep, BinList, etc). + [default: None, meaning ten equally-sized bins] + @param xlabel The label for the x-axis. + [default: None, meaning use the value of `x_field`] + @param ylabel The label for the y-axis. + [default: None, meaning use a string based on the `method` + and `y_value`] + @param color The color of scattered points. This color is also applied to + linear regression if argument `linear_regression` is True. + [default: None, meaning follow a matplotlib's default color] + @param lim The limit of the axes. This can be specified explicitly by + using tuples such as ((xmin, xmax), (ymin, ymax)). + [default: None, meaning do not set any limits] + @param equal_axis If True, force ticks of the x-axis and y-axis equal to each other. + [default: False] + @param linear_regression If True, perform linear regression for x and y and plot a + regression line. If yerr is not None, perform the linear + regression incorporating the error into the standard chi^2 + and plot a regression line with a 1-sigma allowed region. + [default: False] + @param reference_line Draw a reference line. If reference_line == 'one-to-one', x=y + is drawn. If reference_line == 'zero', y=0 is drawn. + A user-specific function can be used by passing an object which + has an attribute '__call__' and returns a 1-d Numpy array. + @returns a matplotlib.figure.Figure object + """ + if not method: + method = self.method + if not x_field: + if self.x_field: + x_field = self.x_field + else: + raise RuntimeError('Must pass x_field kwarg if not defined when initializing object') + if not y_field: + if self.y_field: + y_field = self.y_field + elif not hasattr(method, '__call__') and not method=='count': + raise RuntimeError('Must pass y_field kwarg if not defined when initializing object') + if not yerr_field: + yerr_field = self.yerr_field + if not w_field: + w_field = self.w_field + if w_field and yerr_field: + raise RuntimeError("Cannot pass both a yerr_field and a w_field") + if area is not None and method is not 'count': + raise RuntimeError('Can only pass keyword argument area if method is "count"') + if area_units is not None and area is None: + raise RuntimeError('Can only pass keyword argument area_units if also passing area') + nans = numpy.isnan(array[x_field]) + if y_field: + nans = nans | numpy.isnan(array[y_field]) + if w_field: + nans = nans | numpy.isnan(array[w_field]) + if yerr_field: + nans = nans | numpy.isnan(array[yerr_field]) + print("Skipping", numpy.sum(nans), "nans out of", len(array)) + array = array[numpy.invert(nans)] + if not binning: + binning = self.binning + if not isinstance(binning, + (stile.binning.BinStep, stile.binning.BinList, stile.binning.BinFunction)): + try: + low = min(array[x_field]) + high = max(array[x_field]) + high += 1.E-6*(high-low) # to avoid <= upper bound problem + binning = stile.BinStep(field=x_field, low=low, high=high, n_bins=numpy.round(binning)) + except: + raise RuntimeError("Cannot understand binning argument: %s. Must be a " + "stile.BinStep, stile.BinList, or stile.BinFunction, or " + "a number"%str(binning)) + + x_vals = [] + y_vals = [] + yerr_vals = [] + for ibin, bin in enumerate(binning()): + masked_array = bin(array) + if w_field: + weights = masked_array[w_field] + elif yerr_field: + weights = 1./masked_array[yerr_field]**2 + else: + weights = numpy.ones(masked_array.shape[0]) + x_vals.append(numpy.mean(masked_array[x_field])) + if method=='mean' or method=='rms': + sum_weights = numpy.sum(weights) + mean = numpy.sum(weights*masked_array[y_field])/sum_weights + # This is the unbiased sigma estimator for weighted data, with a further /sqrt(n) + # for error on the mean + sigma = numpy.sqrt(numpy.sum((weights*(masked_array[y_field]-mean))**2)/ + (sum_weights-numpy.sum(weights**2)/sum_weights)) + if method=='mean': + y_vals.append(mean) + yerr_vals.append(sigma/numpy.sqrt(len(masked_array))) + elif method=='rms': + rms = numpy.sqrt(numpy.sum(weights*masked_array[y_field]**2)/ + numpy.sum(weights)) + y_vals.append(rms) + # error estimate for rms: + # sigma_x^2 = - ^2 + # x_{rms}^2 = = sigma_x^2 + ^2 + # x_{rms} = \sqrt{sigma_x^2 + ^2} + # By the Taylor expansion-type propagation of errors, + # sigma_{x_{rms}}^2 = (\partial x_{rms}/\partial sigma_x)^2 sigma_{sigma_x}^2 + + # (\partial x_{rms}/\partial )^2 sigma_{}^2 + # = sigma_x^2/x_rms^2 * sigma_{sigma_x}^2 + ^2/x_rms^2 * sigma_{}^2 + # sigma_{} is sigma_x/sqrt(n), and sigma_{sigma_x} is approx + # 1/sqrt(2(n-1))*sigma_x for large n + # (see http://ai.eecs.umich.edu/~fessler/papers/files/tr/stderr.pdf ). + yerr_vals.append(numpy.sqrt( + (sigma**2/(numpy.sqrt(2*(len(masked_array)-1))*rms))**2 + + (mean*sigma/(numpy.sqrt(len(masked_array))*rms))**2)) + elif method=='median': + y_vals.append(numpy.median(masked_array[y_field])) + # Only correct in the limit of lots of data--but probably where we are + yerr_vals.append( + 1.253*numpy.std(masked_array[y_field])/ + numpy.sqrt(len(masked_array))) + elif method=='count': + # I think there is a better method of weighted counting we could be using? + y_vals.append(numpy.sum(weights)) + yerr_vals.append(numpy.sqrt(numpy.sum(weights))) + if area is not None: + y_vals /= area + yerr_vals /= area + elif hasattr(method, '__call__'): + val = method(masked_array) + if hasattr(val, 'len'): + if len(val)==2: + y_vals.append(val[0]) + yerr_vals.append(val[1]) + else: + raise RuntimeError('method callable may return a tuple of at most length 2 ' + '(a value and its error)') + else: + y_vals.append(val) + if hasattr(method, '__call__'): + y_name = "f(data)" + elif method=='count': + if area: + if area_units: + y_name = "number density [1/%s]"%area_units + else: + y_name = "number density" + else: + y_name = 'number of objects' + else: + y_name = method + " of " + y_field + if isinstance(binning, stile.binning.BinFunction): + x_name = 'bin number' + else: + x_name = x_field + x_vals = numpy.array(x_vals) + y_vals = numpy.array(y_vals) + if yerr_vals: + yerr_vals = numpy.array(yerr_vals) + self.data = numpy.rec.fromarrays([x_vals, y_vals, yerr_vals], + names = [x_name, y_name, y_name+'_err']) + else: + yerr_vals = None + self.data = numpy.rec.fromarrays([x_vals, y_vals], + names = [x_name, y_name]) + if not xlabel: + xlabel = x_name + if not ylabel: + ylabel = y_name + return self.scatterPlot(x_vals, y_vals, yerr_vals, None, + xlabel=xlabel, ylabel=ylabel, + color=color, lim=lim, equal_axis=equal_axis, + linear_regression=linear_regression, reference_line=reference_line) + diff --git a/tests/test_sys_tests.py b/tests/test_sys_tests.py new file mode 100755 index 0000000..3a4d49a --- /dev/null +++ b/tests/test_sys_tests.py @@ -0,0 +1,125 @@ +import numpy +import sys +import unittest + +try: + import stile +except ImportError: + sys.path.append('..') + import stile + + +e1 = [0.421276226426, 0.248889539829, 0.148089741551, 0.566468893815, 0.770334448403, + -1.27622668103, -0.278332559516, -0.0251835420899, -1.62818069301, 0.702161390942, + 0.129529824212, -1.27776942697, -0.710875839398, -1.1335793859, -1.04080707461, + -1.65086191504, -1.40401920141, 0.7444639257, 0.433889348913, -0.41658808983, 0.0595854281922, + -0.954522564497, 0.289112274197, -0.0269082232735, 0.524005376038, -0.884279805103, + 0.386267187801, 0.459108410812, -0.837961549229, -2.94522510662, 0.217256993842, + 1.03461460385, 0.747697610909, -1.04474109671, 0.148850586487, -0.25280413009, 0.439367536488, + -0.654960854765, 0.720404748973, -0.294096536953, 0.743294479812, -0.876907911916, + 1.7968329893, 0.525003247233, 0.843964480056, -0.0439629812406, 0.734188694026, + 0.107428156418, -0.708738286387, 2.10618204709, 0.0375719635484, 1.28394004218, + -0.0148659488638, 0.335083593734, 0.270957283902, -0.0216356036923, 0.102495479916, + 0.0788080513263, -0.152250345965, -1.9232014627, 1.76281941927, -0.77827391566, + -0.327698245345, -2.8349486293, 0.486518651676, -0.174461225907, 1.49712076004, + -0.0192066839702, -0.440034592757, 0.531829410488, 1.00373231941, 1.00209273045, 1.1443289089, + 0.400742577523, 0.339739232257, -1.44558832832, -0.694514979012, -0.745566597959, + 1.75860234203, -0.0884143619186, 1.71907390731, 1.60589446927, 0.30850571315, -2.55379214159, + 0.99217104105, 1.09011037649, 0.0380480602541, -0.724767050203, -0.0229718168813, + -0.123116588353, 1.06084090472, -3.06698994438, 0.753312286232, -1.27553750378, + -1.19888696958, -0.505577200354, -1.53081183926, 0.721787458056, -1.64629351295, + 2.78859472966, -0.628463116903, 1.83803080657, -1.22681682534, 1.53392950221, -0.146953706144, + -0.223805341429, -0.934543873994, -2.47187956559, -0.172048621467, 0.626254208767, + -0.320556115122, -0.152165788827, 0.176846339142, 0.799374722463, -0.528247969287, + 0.0175850586217, 2.29336191268, -1.81488630602, -0.091685410815, 1.94971486187] +mag = [14.8435156914, 14.1908463802, 14.3260978164, 14.7713128658, 14.2375565438, 14.7216183965, + 14.2066642651, 14.6890539665, 14.3540117115, 14.4674828168, 14.2725141398, 14.6878726633, + 14.8357036783, 14.8605710431, 14.2105982155, 14.6843866996, 14.1853052891, 14.7347528304, + 14.260101013, 14.0080025531, 15.5034102914, 15.8914987723, 15.3574767202, 15.2144006756, + 15.3464926008, 15.8491548213, 15.2385574084, 15.32234957, 15.7373386709, 15.3210672383, + 15.3820830471, 15.2247734595, 15.7135322218, 15.5777172989, 15.8207178049, 15.3497158703, + 15.733125029, 15.0961534562, 15.2008993803, 15.3465561104, 16.736433564, 16.9498795809, + 16.5998325489, 16.5615294886, 16.3587625526, 16.4944837469, 16.1678565896, 16.8446540809, + 16.7225771659, 16.2774487571, 16.3307575884, 16.6506168378, 16.2848323207, 16.9597441601, + 16.0599585425, 16.2993441047, 16.5644310218, 16.0404524774, 16.1004008827, 16.105415745, + 17.9515814785, 17.4868439292, 17.7695008586, 17.5840244446, 17.6474927952, 17.1217961689, + 17.1482286314, 17.4023940335, 17.7749584264, 17.6197726047, 17.538029214, 17.1410107284, + 17.2664306017, 17.0677492338, 17.2480659229, 17.8482180569, 17.9289409015, 17.056659366, + 17.8604921386, 17.4703316731, 18.1843444943, 18.5056441016, 18.1254228416, 18.5605439377, + 18.9192087749, 18.9835669315, 18.3701788909, 18.0059356899, 18.5994306164, 18.7364872086, + 18.847013152, 18.2390753809, 18.2296005197, 18.334774316, 18.5773209381, 18.9908561801, + 18.0341188447, 18.0163669231, 18.4686584294, 18.0542893363, 19.4487328726, 19.0490591931, + 19.8500530986, 19.1794731161, 19.562552228, 19.2382319793, 19.2138138099, 19.9348404502, + 19.784382543, 19.2725371274, 19.5729389869, 19.7636863178, 19.3471996974, 19.6713559119, + 19.7096265534, 19.9638882085, 19.1217978937, 19.1296452118, 19.0901038016, 19.3641559185] +data = numpy.rec.fromarrays([numpy.array(mag), numpy.array(e1)], names=['mag','e1']) + +means = numpy.array([-0.33386605345, -0.143461455482, 0.261209398388, + 0.118940939595, -0.119718822072, 0.0503444521599]) +medians = numpy.array([-0.1517580508039, 0.10421800734, 0.104961818167, + 0.160266274143, 0.00753812168641, -0.149559747485]) +rms = numpy.array([0.898951027918, 0.890891207925, 0.911144568394, + 1.11059740575, 1.51276947794, 1.16757299242]) +counts = [20, 20, 20, 20, 18, 22] + +def e1_median(array): + return numpy.median(array['e1']) + +class TestSysTests(unittest.TestCase): + def test_BinnedScatterPlotSysTest(self): + """ + Test that the values passed to scatterPlot from BinnedScatterPlotSysTest.__call__ make sense + """ + test_obj_1 = stile.BinnedScatterPlotSysTest() # Blank, to test kwarg calls + test_obj_2 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6) + test_obj_3 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', + binning=stile.BinStep(field='mag', low=14.0080025531, high=19.9638941644, n_bins=6)) + test_obj_4 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, + method='median') + test_obj_5 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, + method='rms') + test_obj_6 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, + method='count') + test_obj_7 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, + method=e1_median) + + test_obj_1(data, x_field='mag', y_field='e1', binning=6) + mean_obj_1 = test_obj_1.getData() + numpy.testing.assert_almost_equal(mean_obj_1['mean of e1'], means) + test_obj_2(data) + mean_obj_2 = test_obj_2.getData() + numpy.testing.assert_almost_equal(mean_obj_2['mean of e1'], means) + test_obj_3(data) + mean_obj_3 = test_obj_3.getData() + numpy.testing.assert_almost_equal(mean_obj_3['mean of e1'], means) + test_obj_1(data, x_field='mag', y_field='e1', binning=6, method='median') + median_obj_1 = test_obj_1.getData() + numpy.testing.assert_almost_equal(median_obj_1['median of e1'], medians) + test_obj_4(data) + median_obj_4 = test_obj_4.getData() + numpy.testing.assert_almost_equal(median_obj_4['median of e1'], medians) + test_obj_1(data, x_field='mag', y_field='e1', binning=6, method='rms') + rms_obj_1 = test_obj_1.getData() + numpy.testing.assert_almost_equal(rms_obj_1['rms of e1'], rms) + test_obj_5(data) + rms_obj_5 = test_obj_5.getData() + numpy.testing.assert_almost_equal(rms_obj_5['rms of e1'], rms) + + test_obj_1(data, x_field='mag', y_field='e1', binning=6, method='count') + count_obj_1 = test_obj_1.getData() + numpy.testing.assert_equal(count_obj_1['count of e1'], counts) + test_obj_6(data) + count_obj_6 = test_obj_6.getData() + numpy.testing.assert_equal(count_obj_6['count of e1'], counts) + + + test_obj_1(data, x_field='mag', y_field='e1', binning=6, method=e1_median) + median_obj_1 = test_obj_1.getData() + numpy.testing.assert_almost_equal(median_obj_1['f(data)'], medians) + test_obj_7(data) + median_obj_7 = test_obj_7.getData() + numpy.testing.assert_almost_equal(median_obj_7['f(data)'], medians) + +if __name__=='__main__': + unittest.main() + \ No newline at end of file