From ab3bae6592e4d51d9639b44aaa60ab786b214f68 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Sun, 14 Jun 2015 20:31:30 -0400 Subject: [PATCH 01/18] Adding a first run at a BinnedScatterPlot class (#68) --- stile/__init__.py | 2 +- stile/sys_tests.py | 181 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 182 insertions(+), 1 deletion(-) diff --git a/stile/__init__.py b/stile/__init__.py index c32fad1..0b2241a 100644 --- a/stile/__init__.py +++ b/stile/__init__.py @@ -6,4 +6,4 @@ from .treecorr_utils import ReadTreeCorrResultsFile from .data_handler import DataHandler from .sys_tests import (GalaxyShearSysTest, BrightStarShearSysTest, StarXGalaxyShearSysTest, - StarXStarShearSysTest, StatSysTest) + StarXStarShearSysTest, StatSysTest, BinnedScatterPlotSysTest) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index a754d2e..e63e407 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1425,3 +1425,184 @@ def __call__(self, array, per_ccd_stat = None, color = '', lim=None): linear_regression=True, reference_line='zero') + +class BinnedScatterPlotSysTest(ScatterPlotSysTest): + 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. + + 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' or 'rms'), 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 + 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): + """ + 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. + @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' or 'rms'), 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] + @param xlabel The label for the x-axis. + [default: None, meaning do not show a label on the x-axis] + @param ylabel The label for the y-axis. + [default: None, meaning do not show a label on the y-axis] + @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] + @equal_axis If True, force ticks of the x-axis and y-axis equal to each other. + [default: False] + @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] + @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 x_field: + if self.x_field: + x_field = self.x_field + else: + raise ValueError('Must pass x_field kwarg if not defined when initializing object') + if not y_field: + if self.y_field: + y_field = self.y_field + else: + raise ValueError('Must pass x_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 not method: + method = self.method + 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(rray[x_field]) + high += 1.E-6*(high-low) # to avoid <= upper bound problem + binning = stile.BinStep(low=low, high=high, field=x_field, 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)) + if w_field: + if yerr_field: + raise RuntimeError("Cannot pass both a yerr_field and a w_field") + weights = array[w_field] + elif yerr_field: + weights = 1./array[yerr_field]**2 + else: + weights = numpy.ones(array.shape[0]) + x_vals = [] + y_vals = [] + y_err_vals = [] + for ibin, bin in enumerate(binning()): + mask = bin(array) + x_vals.append(numpy.mean(array[x_field][mask])) + if method=='mean': + y_vals.append(numpy.sum(weights*array[y_field])[mask]/numpy.sum(weights[mask])) + yerr_vals.append() # figure this out + elif method=='median': + y_vals.append(numpy.median(array[y_field][mask])) + yerr_vals.append() # figure this out + elif method=='rms': + y_vals.append(numpy.sqrt(numpy.sum((weights*array[y_field]**2)[mask])/numpy.sum(weights[mask]))) + elif hasattr(method, '__call__'): + val = method(array[mask]) + 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)" + 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, y, yerr, None, + xlabel=xlabel, ylabel=ylabel, + color=color, lim=lim, equal_axis=False, + linear_regression=True, reference_line=reference_line) + From 480196a20d06642228886d33b381488f01ef4051 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Sun, 14 Jun 2015 20:32:12 -0400 Subject: [PATCH 02/18] Adding a first run at a testing module for the BinnedScatterPlot (#68) --- tests/test_sys_tests.py | 109 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100755 tests/test_sys_tests.py diff --git a/tests/test_sys_tests.py b/tests/test_sys_tests.py new file mode 100755 index 0000000..b96ddb1 --- /dev/null +++ b/tests/test_sys_tests.py @@ -0,0 +1,109 @@ +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.329512262063, -0.143461455482, 0.261209398388, + 0.118940939595, -0.119718822072, 0.0503444521599]) +medians = numpy.array([-0.0251835420899, 0.10421800734, 0.104961818167, + 0.160266274143, 0.00753812168641, -0.149559747485]) +rms = numpy.array([0.917339212222, 0.890891207925, 0.911144568394, + 1.11059740575, 1.51276947794, 1.16757299242]) + +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(low=14.0080025531, high=9.9638941644, n_bins=6)) + test_obj_4 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, method='median') + test_obj_4 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, method='rms') + test_obj_4 = 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_equal(mean_obj_1['e1'], means) + test_obj_2(data) + mean_obj_2 = test_obj_2.getData() + numpy.testing.assert_equal(mean_obj_2['e1'], means) + test_obj_3(data) + mean_obj_3 = test_obj_3.getData() + numpy.testing.assert_equal(mean_obj_3['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_equal(median_obj_1['e1'], medians) + test_obj_4(data) + median_obj_4 = test_obj_4.getData() + numpy.testing.assert_equal(median_obj_4['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_equal(rms_obj_1['e1'], rms) + test_obj_5(data) + rms_obj_5 = test_obj_5.getData() + numpy.testing.assert_equal(rms_obj_5['e1'], rms) + 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_equal(median_obj_1['e1'], medians) + test_obj_6(data) + median_obj_6 = test_obj_6.getData() + numpy.testing.assert_equal(median_obj_6['e1'], medians) + +if __name__=='__main__': + unittest.main() + \ No newline at end of file From dd61de182b13ec2609c7500327710bcb4447aa5a Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Tue, 16 Jun 2015 01:11:54 -0400 Subject: [PATCH 03/18] Bugfixes found in running tests (#68) --- stile/sys_tests.py | 50 +++++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index e63e407..d7cbe9b 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1537,39 +1537,39 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu binning = self.binning if not isinstance(binning, (stile.binning.BinStep, stile.binning.BinList, stile.binning.BinFunction)): - try: +# try: low = min(array[x_field]) - high = max(rray[x_field]) + high = max(array[x_field]) high += 1.E-6*(high-low) # to avoid <= upper bound problem - binning = stile.BinStep(low=low, high=high, field=x_field, 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)) - if w_field: - if yerr_field: - raise RuntimeError("Cannot pass both a yerr_field and a w_field") - weights = array[w_field] - elif yerr_field: - weights = 1./array[yerr_field]**2 - else: - weights = numpy.ones(array.shape[0]) + 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)) + if w_field and yerr_field: + raise RuntimeError("Cannot pass both a yerr_field and a w_field") x_vals = [] y_vals = [] - y_err_vals = [] + yerr_vals = [] for ibin, bin in enumerate(binning()): - mask = bin(array) - x_vals.append(numpy.mean(array[x_field][mask])) + 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': - y_vals.append(numpy.sum(weights*array[y_field])[mask]/numpy.sum(weights[mask])) - yerr_vals.append() # figure this out + y_vals.append(numpy.sum(weights*masked_array[y_field])/numpy.sum(weights)) + #yerr_vals.append() # figure this out elif method=='median': - y_vals.append(numpy.median(array[y_field][mask])) - yerr_vals.append() # figure this out + y_vals.append(numpy.median(masked_array[y_field])) + #yerr_vals.append() # figure this out elif method=='rms': - y_vals.append(numpy.sqrt(numpy.sum((weights*array[y_field]**2)[mask])/numpy.sum(weights[mask]))) + y_vals.append(numpy.sqrt(numpy.sum(weights*masked_array[y_field]**2)/numpy.sum(weights))) elif hasattr(method, '__call__'): - val = method(array[mask]) + val = method(masked_array) if hasattr(val, 'len'): if len(val)==2: y_vals.append(val[0]) @@ -1601,7 +1601,7 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu xlabel = x_name if not ylabel: ylabel = y_name - return self.scatterPlot(x, y, yerr, None, + return self.scatterPlot(x_vals, y_vals, yerr_vals, None, xlabel=xlabel, ylabel=ylabel, color=color, lim=lim, equal_axis=False, linear_regression=True, reference_line=reference_line) From 578430e280e63209e75ba3a71cdf783d898b40d1 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Tue, 16 Jun 2015 01:12:15 -0400 Subject: [PATCH 04/18] Bugfixes to test_sys_tests.py (#68)y --- tests/test_sys_tests.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/tests/test_sys_tests.py b/tests/test_sys_tests.py index b96ddb1..3bb0d78 100755 --- a/tests/test_sys_tests.py +++ b/tests/test_sys_tests.py @@ -54,13 +54,13 @@ 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.329512262063, -0.143461455482, 0.261209398388, +means = numpy.array([-0.33386605345, -0.143461455482, 0.261209398388, 0.118940939595, -0.119718822072, 0.0503444521599]) -medians = numpy.array([-0.0251835420899, 0.10421800734, 0.104961818167, +medians = numpy.array([-0.1517580508039, 0.10421800734, 0.104961818167, 0.160266274143, 0.00753812168641, -0.149559747485]) -rms = numpy.array([0.917339212222, 0.890891207925, 0.911144568394, +rms = numpy.array([0.898951027918, 0.890891207925, 0.911144568394, 1.11059740575, 1.51276947794, 1.16757299242]) - + def e1_median(array): return numpy.median(array['e1']) @@ -71,38 +71,38 @@ def test_BinnedScatterPlotSysTest(self): """ 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(low=14.0080025531, high=9.9638941644, n_bins=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_4 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, method='rms') - test_obj_4 = stile.BinnedScatterPlotSysTest(x_field='mag', y_field='e1', binning = 6, method=e1_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=e1_median) test_obj_1(data, x_field='mag', y_field='e1', binning=6) mean_obj_1 = test_obj_1.getData() - numpy.testing.assert_equal(mean_obj_1['e1'], means) + 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_equal(mean_obj_2['e1'], means) + 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_equal(mean_obj_3['e1'], means) + 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_equal(median_obj_1['e1'], medians) + 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_equal(median_obj_4['e1'], medians) + 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_equal(rms_obj_1['e1'], rms) + 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_equal(rms_obj_5['e1'], rms) + 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=e1_median) median_obj_1 = test_obj_1.getData() - numpy.testing.assert_equal(median_obj_1['e1'], medians) + numpy.testing.assert_almost_equal(median_obj_1['f(data)'], medians) test_obj_6(data) median_obj_6 = test_obj_6.getData() - numpy.testing.assert_equal(median_obj_6['e1'], medians) + numpy.testing.assert_almost_equal(median_obj_6['f(data)'], medians) if __name__=='__main__': unittest.main() From 65674994f92bde41567dcfcc8777d2f381bbc174 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Wed, 17 Jun 2015 18:58:45 -0400 Subject: [PATCH 05/18] Add error calculations and a "count" option to BinnedScatterPlot (#68) --- stile/sys_tests.py | 53 +++++++++++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index d7cbe9b..e67e121 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1446,10 +1446,10 @@ class BinnedScatterPlotSysTest(ScatterPlotSysTest): @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' or 'rms'), 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'] + 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 @@ -1484,10 +1484,10 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu @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' or 'rms'), 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'] + 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 @@ -1537,15 +1537,15 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu binning = self.binning if not isinstance(binning, (stile.binning.BinStep, stile.binning.BinList, stile.binning.BinFunction)): -# try: + 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)) + except: + raise RuntimeError("Cannot understand binning argument: %s. Must be a " + "stile.BinStep, stile.BinList, or stile.BinFunction, or " + "a number"%str(binning)) if w_field and yerr_field: raise RuntimeError("Cannot pass both a yerr_field and a w_field") x_vals = [] @@ -1561,13 +1561,28 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu weights = numpy.ones(masked_array.shape[0]) x_vals.append(numpy.mean(masked_array[x_field])) if method=='mean': - y_vals.append(numpy.sum(weights*masked_array[y_field])/numpy.sum(weights)) - #yerr_vals.append() # figure this out + sum_weights = numpy.sum(weights) + mean = numpy.sum(weights*masked_array[y_field])/sum_weights + y_vals.append(mean) + # This is the unbiased sigma estimator for weighted data, with a further /sqrt(n) + # for error on the mean + yerr_vals.append( + numpy.sqrt(numpy.sum((weights*(masked_array[y_field]-mean))**2)/ + (sum_weights-numpy.sum(weights**2)/sum_weights)/len(masked_array))) elif method=='median': y_vals.append(numpy.median(masked_array[y_field])) - #yerr_vals.append() # figure this out + # 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=='rms': - y_vals.append(numpy.sqrt(numpy.sum(weights*masked_array[y_field]**2)/numpy.sum(weights))) + y_vals.append( + numpy.sqrt(numpy.sum(weights*masked_array[y_field]**2)/ + numpy.sum(weights))) + elif method=='count': + # I think there is a better method of weight-counting we could be using? + y_vals.append(numpy.sum(weights)) + yerr_vals.append(numpy.sqrt(numpy.sum(weights))) elif hasattr(method, '__call__'): val = method(masked_array) if hasattr(val, 'len'): @@ -1603,6 +1618,6 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu ylabel = y_name return self.scatterPlot(x_vals, y_vals, yerr_vals, None, xlabel=xlabel, ylabel=ylabel, - color=color, lim=lim, equal_axis=False, - linear_regression=True, reference_line=reference_line) + color=color, lim=lim, equal_axis=equal_axis, + linear_regression=linear_regression, reference_line=reference_line) From 28c71bba76bc01a1ea473defc3d69056bb413404 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Wed, 17 Jun 2015 18:59:46 -0400 Subject: [PATCH 06/18] Edit line lengths in test_sys_tests.py (#68) --- tests/test_sys_tests.py | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/tests/test_sys_tests.py b/tests/test_sys_tests.py index 3bb0d78..3a4d49a 100755 --- a/tests/test_sys_tests.py +++ b/tests/test_sys_tests.py @@ -60,6 +60,7 @@ 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']) @@ -71,10 +72,16 @@ def test_BinnedScatterPlotSysTest(self): """ 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=e1_median) + 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() @@ -97,12 +104,21 @@ def test_BinnedScatterPlotSysTest(self): 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_6(data) - median_obj_6 = test_obj_6.getData() - numpy.testing.assert_almost_equal(median_obj_6['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() From 518e2ff0bd62c02843728132aeeed2a6a80b67fb Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Wed, 17 Jun 2015 22:26:51 -0400 Subject: [PATCH 07/18] Implement RMS error + add more documentation (#68) --- stile/sys_tests.py | 58 ++++++++++++++++++++++++++++++++++++---------- 1 file changed, 46 insertions(+), 12 deletions(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index e67e121..a5ec5a1 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1436,6 +1436,20 @@ class BinnedScatterPlotSysTest(ScatterPlotSysTest): 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 @@ -1475,6 +1489,10 @@ def __call__(self, array, x_field=None, y_field=None, yerr_field=None, w_field=N 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 @@ -1525,8 +1543,8 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu if not y_field: if self.y_field: y_field = self.y_field - else: - raise ValueError('Must pass x_field kwarg if not defined when initializing object') + elif not hasattr(method, '__call__'): + raise ValueError('Must pass y_field kwarg if not defined when initializing object') if not yerr_field: yerr_field = self.yerr_field if not w_field: @@ -1560,27 +1578,43 @@ class and return the `matplotlib.figure.Figure` that this `__call__` method retu else: weights = numpy.ones(masked_array.shape[0]) x_vals.append(numpy.mean(masked_array[x_field])) - if method=='mean': + if method=='mean' or method=='rms': sum_weights = numpy.sum(weights) mean = numpy.sum(weights*masked_array[y_field])/sum_weights - y_vals.append(mean) # This is the unbiased sigma estimator for weighted data, with a further /sqrt(n) # for error on the mean - yerr_vals.append( - numpy.sqrt(numpy.sum((weights*(masked_array[y_field]-mean))**2)/ - (sum_weights-numpy.sum(weights**2)/sum_weights)/len(masked_array))) + 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=='rms': - y_vals.append( - numpy.sqrt(numpy.sum(weights*masked_array[y_field]**2)/ - numpy.sum(weights))) + elif method=='count': - # I think there is a better method of weight-counting we could be using? + # 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))) elif hasattr(method, '__call__'): From 1834c7c40839248ee340d79e1858d67979aae97e Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Wed, 17 Jun 2015 22:59:16 -0400 Subject: [PATCH 08/18] Add HSC adapters for BinnedScatterPlot specific instances (#68) --- stile/hsc/base_tasks.py | 7 +++--- stile/hsc/sys_test_adapters.py | 46 ++++++++++++++++++++++++++++++++++ stile/sys_tests.py | 5 ++++ 3 files changed, 55 insertions(+), 3 deletions(-) diff --git a/stile/hsc/base_tasks.py b/stile/hsc/base_tasks.py index 04f44c9..713796c 100644 --- a/stile/hsc/base_tasks.py +++ b/stile/hsc/base_tasks.py @@ -63,7 +63,8 @@ class CCDSingleEpochStileConfig(lsst.pex.config.Config): "WhiskerPlotStar", "WhiskerPlotPSF", "WhiskerPlotResidual", "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", - "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma" + "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", + "RMSE1Sky", "RMSE2Sky," "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control TreeCorr", keytype=str, itemtype=str, @@ -231,12 +232,12 @@ def removeFlaggedObjects(self, catalog): @returns The source catalog, masked to the rows which don't have any of our defined flags set. """ - masks = [] + masks = [] if self.config.flags_keep_false: masks += [catalog[flag]==False for flag in self.config.flags_keep_false] if self.config.flags_keep_true: masks += [catalog[flag]==True for flag in self.config.flags_keep_true] - if masks: + if masks: mask = masks[0] for new_mask in masks[1:]: mask = numpy.logical_and(mask, new_mask) diff --git a/stile/hsc/sys_test_adapters.py b/stile/hsc/sys_test_adapters.py index 70b7826..b76e4ea 100644 --- a/stile/hsc/sys_test_adapters.py +++ b/stile/hsc/sys_test_adapters.py @@ -459,6 +459,47 @@ 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 RMSE1Sky(ShapeSysTestAdapter): + 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 RMSE2Sky(ShapeSysTestAdapter): + 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 RMSE1Chip(ShapeSysTestAdapter): + 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 RMSE2Chip(ShapeSysTestAdapter): + 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 CountPerMagnitude(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']) + + + adapter_registry.register("StatsPSFFlux", StatsPSFFluxAdapter) adapter_registry.register("GalaxyShear", GalaxyShearAdapter) adapter_registry.register("BrightStarShear", BrightStarShearAdapter) @@ -474,3 +515,8 @@ def __call__(self,task_config,*data, **kwargs): adapter_registry.register("ScatterPlotResidualVsPSFG1", ScatterPlotResidualVsPSFG1Adapter) adapter_registry.register("ScatterPlotResidualVsPSFG2", ScatterPlotResidualVsPSFG2Adapter) adapter_registry.register("ScatterPlotResidualVsPSFSigma", ScatterPlotResidualVsPSFSigmaAdapter) +adapter_registry.register("RMSE1Sky", RMSE1SkyAdapter) +adapter_registry.register("RMSE2Sky", RMSE2SkyAdapter) +adapter_registry.register("RMSE1Chip", RMSE1ChipAdapter) +adapter_registry.register("RMSE2Chip", RMSE2ChipAdapter) +adapter_registry.register("CountPerMagnitude", CountPerMagnitude) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index a5ec5a1..ccdfd1c 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1477,6 +1477,11 @@ def __init__(self, x_field=None, y_field=None, yerr_field=None, w_field=None, me 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 self.method = method self.binning = binning From 91401362be89ddab6db82504e7f9282de94772b1 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Thu, 2 Jul 2015 03:15:45 -0400 Subject: [PATCH 09/18] Bugfixes making #68 work --- stile/hsc/base_tasks.py | 11 +++++++---- stile/hsc/sys_test_adapters.py | 12 ++++++------ stile/sys_tests.py | 3 ++- 3 files changed, 15 insertions(+), 11 deletions(-) diff --git a/stile/hsc/base_tasks.py b/stile/hsc/base_tasks.py index 713796c..818652b 100644 --- a/stile/hsc/base_tasks.py +++ b/stile/hsc/base_tasks.py @@ -64,7 +64,7 @@ class CCDSingleEpochStileConfig(lsst.pex.config.Config): "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", - "RMSE1Sky", "RMSE2Sky," "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" + "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control TreeCorr", keytype=str, itemtype=str, @@ -658,7 +658,8 @@ class VisitSingleEpochStileConfig(CCDSingleEpochStileConfig): "WhiskerPlotStar", "WhiskerPlotPSF", "WhiskerPlotResidual", "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", - "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma" + "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", + "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control treecorr", keytype=str, itemtype=str, @@ -870,7 +871,8 @@ class PatchSingleEpochStileConfig(CCDSingleEpochStileConfig): "WhiskerPlotStar", "WhiskerPlotPSF", "WhiskerPlotResidual", "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", - "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma" + "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", + "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" ]) whiskerplot_figsize = lsst.pex.config.ListField(dtype=float, doc="figure size for whisker plot", default = [7., 5.]) @@ -953,7 +955,8 @@ class TractSingleEpochStileConfig(CCDSingleEpochStileConfig): "WhiskerPlotStar", "WhiskerPlotPSF", "WhiskerPlotResidual", "ScatterPlotStarVsPSFG1", "ScatterPlotStarVsPSFG2", "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", - "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma" + "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", + "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control treecorr", keytype=str, itemtype=str, diff --git a/stile/hsc/sys_test_adapters.py b/stile/hsc/sys_test_adapters.py index b76e4ea..4035b81 100644 --- a/stile/hsc/sys_test_adapters.py +++ b/stile/hsc/sys_test_adapters.py @@ -459,7 +459,7 @@ 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 RMSE1Sky(ShapeSysTestAdapter): +class RMSE1SkyAdapter(ShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'sky' self.config = config @@ -467,7 +467,7 @@ def __init__(self, config): self.name = 'rms_e1_sky' self.setupMasks(['galaxy']) -class RMSE2Sky(ShapeSysTestAdapter): +class RMSE2SkyAdapter(ShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'sky' self.config = config @@ -475,7 +475,7 @@ def __init__(self, config): self.name = 'rms_e2_sky' self.setupMasks(['galaxy']) -class RMSE1Chip(ShapeSysTestAdapter): +class RMSE1ChipAdapter(ShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'chip' self.config = config @@ -483,7 +483,7 @@ def __init__(self, config): self.name = 'rms_e1_chip' self.setupMasks(['galaxy']) -class RMSE2Chip(ShapeSysTestAdapter): +class RMSE2ChipAdapter(ShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'chip' self.config = config @@ -491,7 +491,7 @@ def __init__(self, config): self.name = 'rms_e2_chip' self.setupMasks(['galaxy']) -class CountPerMagnitude(BaseSysTestAdapter): +class CountPerMagnitudeAdapter(BaseSysTestAdapter): def __init__(self, config): self.config = config self.sys_test = sys_tests.BinnedScatterPlotSysTest(x_field='mag', w_field='w', method='count') @@ -519,4 +519,4 @@ def __init__(self, config): adapter_registry.register("RMSE2Sky", RMSE2SkyAdapter) adapter_registry.register("RMSE1Chip", RMSE1ChipAdapter) adapter_registry.register("RMSE2Chip", RMSE2ChipAdapter) -adapter_registry.register("CountPerMagnitude", CountPerMagnitude) +adapter_registry.register("CountPerMagnitude", CountPerMagnitudeAdapter) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index ccdfd1c..ec7595c 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1481,7 +1481,8 @@ def __init__(self, x_field=None, y_field=None, yerr_field=None, w_field=None, me for quantity in [x_field, y_field, yerr_field, w_field]: if quantity: list_of_quantities.append(quantity) - self.required_quantities = list_of_quantities + self.required_quantities = [list_of_quantities] # need a list per object type + print self.required_quantities self.method = method self.binning = binning From 8541c76c194019365b218ce2fc2c84f12614dc1f Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Wed, 8 Jul 2015 13:01:03 -0400 Subject: [PATCH 10/18] Fix accidental passing of config in shape adapters (#68) --- stile/hsc/sys_test_adapters.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/stile/hsc/sys_test_adapters.py b/stile/hsc/sys_test_adapters.py index 4035b81..02b0d33 100644 --- a/stile/hsc/sys_test_adapters.py +++ b/stile/hsc/sys_test_adapters.py @@ -186,7 +186,7 @@ def __call__(self, task_config, *data, **kwargs): sys_test itself returns. """ new_data = [self.fixArray(d) for d in data] - return self.sys_test(task_config.corr2_kwargs, *new_data, **kwargs) + return self.sys_test(config=task_config.treecorr_kwargs, *new_data, **kwargs) class GalaxyShearAdapter(ShapeSysTestAdapter): """ @@ -459,7 +459,19 @@ 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(ShapeSysTestAdapter): +class NoConfigShapeSysTestAdapter(ShapeSysTestAdapter): + """ + Just like the ShapeSysTestAdapter, only don't send the config object to the test. + """ + def __call__(self, task_config, *data, **kwargs): + """ + Call this object's sys_test with the given data and kwargs, and return whatever the + sys_test itself returns. + """ + new_data = [self.fixArray(d) for d in data] + return self.sys_test(*new_data, **kwargs) + +class RMSE1SkyAdapter(NoConfigShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'sky' self.config = config @@ -467,7 +479,7 @@ def __init__(self, config): self.name = 'rms_e1_sky' self.setupMasks(['galaxy']) -class RMSE2SkyAdapter(ShapeSysTestAdapter): +class RMSE2SkyAdapter(NoConfigShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'sky' self.config = config @@ -475,7 +487,7 @@ def __init__(self, config): self.name = 'rms_e2_sky' self.setupMasks(['galaxy']) -class RMSE1ChipAdapter(ShapeSysTestAdapter): +class RMSE1ChipAdapter(NoConfigShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'chip' self.config = config @@ -483,7 +495,7 @@ def __init__(self, config): self.name = 'rms_e1_chip' self.setupMasks(['galaxy']) -class RMSE2ChipAdapter(ShapeSysTestAdapter): +class RMSE2ChipAdapter(NoConfigShapeSysTestAdapter): def __init__(self, config): self.shape_type = 'chip' self.config = config From 713b34f4c02b308b2f4a0b1f1a994c186a1e5d24 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Wed, 8 Jul 2015 13:01:45 -0400 Subject: [PATCH 11/18] Fix implementation of count method, plus checking for NaNs (#68) --- stile/sys_tests.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index ec7595c..6e15163 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1541,6 +1541,8 @@ class function `scatterPlot()` unaltered. 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 @@ -1549,14 +1551,12 @@ class function `scatterPlot()` unaltered. if not y_field: if self.y_field: y_field = self.y_field - elif not hasattr(method, '__call__'): + elif not hasattr(method, '__call__') and not method=='count': raise ValueError('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 not method: - method = self.method if not binning: binning = self.binning if not isinstance(binning, @@ -1575,6 +1575,11 @@ class function `scatterPlot()` unaltered. x_vals = [] y_vals = [] yerr_vals = [] + if y_field: + nans = numpy.isnan(array[y_field]) + print "Skipping", numpy.sum(nans), "nans out of", len(array) + array = array[numpy.invert(nans)] + for ibin, bin in enumerate(binning()): masked_array = bin(array) if w_field: @@ -1636,6 +1641,8 @@ class function `scatterPlot()` unaltered. y_vals.append(val) if hasattr(method, '__call__'): y_name = "f(data)" + elif method=='count': + y_name = 'number of objects' else: y_name = method + " of " + y_field if isinstance(binning, stile.binning.BinFunction): From 21248188eae20637bcafb2a61c6ad185776aaae8 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Wed, 8 Jul 2015 13:02:03 -0400 Subject: [PATCH 12/18] Move to getting key first for magnitude stuff (#68) --- stile/hsc/base_tasks.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stile/hsc/base_tasks.py b/stile/hsc/base_tasks.py index 818652b..909ec4a 100644 --- a/stile/hsc/base_tasks.py +++ b/stile/hsc/base_tasks.py @@ -582,7 +582,8 @@ def computeExtraColumn(self, col, data, calib_data, calib_type, xy0=None): elif calib_type=="calexp": zeropoint = 2.5*numpy.log10(calib_data.get("FLUXMAG0")) key = data.schema.find('flux.psf.flags').key - return (zeropoint - 2.5*numpy.log10(data.getPsfFlux()), + flux_key = data.schema.find('flux.psf').key + return (zeropoint - 2.5*numpy.log10([src.get(flux_key) for src in data]), numpy.array([src.get(key)==0 for src in data])) elif col=="w": # Use uniform weights for now if we don't use shapes ("w" will be removed from the From 15cb43d7b9b8be42288f480f142f6fa55dde612f Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Thu, 16 Jul 2015 00:35:14 -0400 Subject: [PATCH 13/18] Check for NaNs in all fields (#68) --- stile/sys_tests.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index 6e15163..3bbc7db 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1482,7 +1482,6 @@ def __init__(self, x_field=None, y_field=None, yerr_field=None, w_field=None, me if quantity: list_of_quantities.append(quantity) self.required_quantities = [list_of_quantities] # need a list per object type - print self.required_quantities self.method = method self.binning = binning @@ -1557,6 +1556,17 @@ class function `scatterPlot()` unaltered. 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") + 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, @@ -1570,16 +1580,10 @@ class function `scatterPlot()` unaltered. raise RuntimeError("Cannot understand binning argument: %s. Must be a " "stile.BinStep, stile.BinList, or stile.BinFunction, or " "a number"%str(binning)) - if w_field and yerr_field: - raise RuntimeError("Cannot pass both a yerr_field and a w_field") + x_vals = [] y_vals = [] yerr_vals = [] - if y_field: - nans = numpy.isnan(array[y_field]) - print "Skipping", numpy.sum(nans), "nans out of", len(array) - array = array[numpy.invert(nans)] - for ibin, bin in enumerate(binning()): masked_array = bin(array) if w_field: From 85cfb2ba048c359dc436d8f778c568a9fc2d01bb Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Thu, 16 Jul 2015 00:43:18 -0400 Subject: [PATCH 14/18] Documentation updates (#68) --- stile/sys_tests.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index 3bbc7db..c333e3a 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1088,14 +1088,14 @@ def scatterPlot(self, x, y, yerr=None, z=None, xlabel=None, ylabel=None, zlabel= If one passes float p, it calculate p%-percentile around median for each of x-axis and y-axis. [default: None, meaning do not set any limits] - @equal_axis If True, force ticks of x-axis and y-axis equal to each other. + @param equal_axis If True, force ticks of x-axis and y-axis equal to each other. [default: False] - @linear_regression If True, perform linear regression for x and y and plot a regression + @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 with incorporating the error into the standard chi^2 and plot a regression line with a 1-sigma allowed region. [default: False] - @reference_line Draw a reference line. If reference_line == 'one-to-one', x=y is + @param reference_line Draw a reference line. If reference_line == 'one-to-one', x=y is drawn. If reference_line == 'zero', y=0 id drawn. A user-specific function can be used by passing an object which has an attribute '__call__' and returns a 1-d Numpy array. @@ -1518,23 +1518,24 @@ class function `scatterPlot()` unaltered. Stile (BinStep, BinList, etc). [default: None, meaning ten equally-sized bins] @param xlabel The label for the x-axis. - [default: None, meaning do not show a label on the x-axis] + [default: None, meaning use the value of `x_field`] @param ylabel The label for the y-axis. - [default: None, meaning do not show a label on 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] - @equal_axis If True, force ticks of the x-axis and y-axis equal to each other. + @param equal_axis If True, force ticks of the x-axis and y-axis equal to each other. [default: False] - @linear_regression If True, perform linear regression for x and y and plot a + @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] - @reference_line Draw a reference line. If reference_line == 'one-to-one', x=y + @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. From 320f3e3a9c9b9a93d8aa1b5a9883a444b8b223ac Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Thu, 5 Nov 2015 06:36:30 -0500 Subject: [PATCH 15/18] Add area kwargs for number density option to BinnedScatterPlot (#68 --- stile/sys_tests.py | 38 +++++++++++++++++++++++++++----------- 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index 2eb3b37..0b9700a 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -1543,7 +1543,8 @@ def __init__(self, x_field=None, y_field=None, yerr_field=None, w_field=None, me 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): + 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 @@ -1566,7 +1567,10 @@ class function `scatterPlot()` unaltered. 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'] + 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 @@ -1603,27 +1607,31 @@ class function `scatterPlot()` unaltered. if self.x_field: x_field = self.x_field else: - raise ValueError('Must pass x_field kwarg if not defined when initializing object') + 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 ValueError('Must pass y_field kwarg if not defined when initializing object') + 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 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)] + array = array[numpy.invert(nans)] if not binning: binning = self.binning if not isinstance(binning, @@ -1637,7 +1645,7 @@ class function `scatterPlot()` unaltered. 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 = [] @@ -1684,11 +1692,13 @@ class function `scatterPlot()` unaltered. 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'): @@ -1702,8 +1712,14 @@ class function `scatterPlot()` unaltered. y_vals.append(val) if hasattr(method, '__call__'): y_name = "f(data)" - elif method=='count': - y_name = 'number of objects' + 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): From b187d741a506b9d9240927155805b05fdf964aa4 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Thu, 5 Nov 2015 17:13:48 -0500 Subject: [PATCH 16/18] Add RMSESky, RMSEChip, and change default tests in hsc module (#68) --- stile/hsc/base_tasks.py | 8 ++++---- stile/hsc/sys_test_adapters.py | 36 +++++++++++++++++++++++++++++++++- 2 files changed, 39 insertions(+), 5 deletions(-) diff --git a/stile/hsc/base_tasks.py b/stile/hsc/base_tasks.py index 42bfcb1..ff84ffd 100644 --- a/stile/hsc/base_tasks.py +++ b/stile/hsc/base_tasks.py @@ -74,7 +74,7 @@ class CCDSingleEpochStileConfig(lsst.pex.config.Config): "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", "ScatterPlotResidualSigmaVsPSFMag", - "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" + "RMSESky", "RMSEChip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control TreeCorr", keytype=str, itemtype=str, @@ -744,7 +744,7 @@ class VisitSingleEpochStileConfig(CCDSingleEpochStileConfig): "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", "ScatterPlotResidualSigmaVsPSFMag", - "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" + "RMSESky", "RMSEChip", "CountPerMagnitude" ]) treecorr_kwargs = lsst.pex.config.DictField(doc="extra kwargs to control treecorr", keytype=str, itemtype=str, @@ -1059,7 +1059,7 @@ class PatchSingleEpochStileConfig(CCDSingleEpochStileConfig): "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", "ScatterPlotResidualSigmaVsPSFMag", - "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" + "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", @@ -1156,7 +1156,7 @@ class TractSingleEpochStileConfig(CCDSingleEpochStileConfig): "ScatterPlotStarVsPSFSigma", "ScatterPlotResidualVsPSFG1", "ScatterPlotResidualVsPSFG2", "ScatterPlotResidualVsPSFSigma", "ScatterPlotResidualSigmaVsPSFMag", - "RMSE1Sky", "RMSE2Sky", "RMSE1Chip", "RMSE2Chip", "CountPerMagnitude" + "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 1e66dd8..f3efbcb 100644 --- a/stile/hsc/sys_test_adapters.py +++ b/stile/hsc/sys_test_adapters.py @@ -569,7 +569,39 @@ def __init__(self, config): 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) @@ -594,3 +626,5 @@ def __init__(self, config): 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) From d7ef8e8abd84485941583a12cc21c090f49addc3 Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Thu, 12 Jul 2018 14:09:57 -0700 Subject: [PATCH 17/18] Fix bug in BinnedScatterPlot inheritance (#68) --- CHANGELOG.md | 1 + stile/sys_tests.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) 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/sys_tests.py b/stile/sys_tests.py index d9e033b..172b6be 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -2362,7 +2362,7 @@ 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(ScatterPlotSysTest): +class BinnedScatterPlotSysTest(BaseScatterPlotSysTest): short_name = 'binnedscatterplot' """ A base class for Stile systematics tests that generate scatter plots binned in one variable, From fb0eaacac0461763e418eb2b65ddd8427ce31a7c Mon Sep 17 00:00:00 2001 From: Melanie Simet Date: Thu, 12 Jul 2018 14:17:32 -0700 Subject: [PATCH 18/18] Fix Python 3 printing bug in BinnedScatterPlotSysTest (#68) --- stile/sys_tests.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index 172b6be..10d80d9 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -2510,7 +2510,7 @@ class function `scatterPlot()` unaltered. 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) + print("Skipping", numpy.sum(nans), "nans out of", len(array)) array = array[numpy.invert(nans)] if not binning: binning = self.binning