diff --git a/CHANGELOG.md b/CHANGELOG.md index 0ac897b..70c2db8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,4 @@ 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) 2/10/15: Add a setup.py installer (issue #57) 8/26/14: Change from relying on compiled C-code corr2 to Python package TreeCorr (issue #33) diff --git a/stile/sys_tests.py b/stile/sys_tests.py index 0dbe9c0..077306b 100644 --- a/stile/sys_tests.py +++ b/stile/sys_tests.py @@ -193,6 +193,33 @@ def CorrelationFunctionSysTest(type=None): type given as the first argument) and that it lacks many of the convenience variables the other CorrelationFunctions have, such as self.objects_list and self.required_quantities. + These produce different estimators depending on the type. Point-point estimates by default use + the Landy-Szalay estimator: + xi = (DD-2DR+RR)/RR + and must include a random catalog. + + All shears in the following descriptions are the complex form in the frame aligned with the + vector between the two points, that is, g = gamma_t + i*gamma_x. + + Point-shear estimates are equivalent to average tangential shear, returned as real + and imaginary . If random catalogs are given, they are used as random lenses, and + data*shear-random*shear is returned instead. + + Shear-shear correlation functions are xi_+ and xi_-. Xi_+ is, nominally, g1 times g2*, and is + given as both the real and imaginary components. xi_+,im should be consistent with 0 to within + noise. Xi_- on the other hand is g1 times g2 (not complex conjugate). Similarly, xi_-,im should + be 0. (Note that g1 and g2 here are two *complex* shears g1_t + i*g1_x and g1_t+i*g2_x from the + two catalogs, not the two components of a single shear in sky coordinates or chip frame.) + + Point-scalar (point-kappa) estimates are equivalent to ; scalar-shear estimates are + equivalent to with a corresponding imaginary case; scalar-scalar estimates + are . Random catalogs result in compensated estimators as in the point-shear + case. + + Aperture mass statistics of various kinds are also available via the + BaseCorrelationFunctionSysTest class; as we do not implement those for any standard Stile tests, + interested users are directed to the TreeCorr documentation for further information. + Users who are writing their own code can of course pass other data types to these functions than the ones given (eg sending galaxy data to a StarXStarSize correlation function); we include options for A) automatic processing, B) ease of understanding code, and C) suggesting @@ -281,7 +308,7 @@ class method :func:`getCF()` which runs a TreeCorr correlation function on a giv PlotDetails(t_field='', t_title=r'$\langle \kappa \rangle$', datarandom_t_field='kappa_', datarandom_t_title='$kappa_{', sigma_field='sigma', y_title="$\kappa$"), # nk - PlotDetails(t_field='xi', t_title=r'$\xi$', sigma_field='sigma_xi', y_title=r"$\xi$"), # k2 + PlotDetails(t_field='xi', t_title=r'$\xi_{\mathrm{re}}$', sigma_field='sigma_xi', y_title=r"$\xi$"), # k2 PlotDetails(t_field='', t_title=r'$\langle \kappa \gamma_T\rangle$', x_field='', x_title=r'$\langle \kappa \gamma_X\rangle$', datarandom_t_field='kgamT_', datarandom_t_title=r'$\kappa \gamma_{T', @@ -516,6 +543,10 @@ def getCF(self, correlation_function_type, data, data2=None, results = stile.ReadTreeCorrResultsFile(output_file) os.close(handle) os.remove(output_file) + names = results.dtype.names + # Add the sep units to the column names of radial bins from TreeCorr outputs + names = [n+' [%s]'%treecorr_kwargs['sep_units'] if 'R' in n else n for n in names] + results.dtype.names = names return results def compensateDefault(self, data, data2, random, random2, both=False): @@ -596,6 +627,7 @@ def plot(self, data, colors=['r', 'b'], log_yscale=False, # Plot the first thing curr_plot = 0 ax = fig.add_subplot(nrows, 1, 1) + ax.axhline(0, alpha=0.7, color='gray') ax.errorbar(data[r], data[pd.t_field], yerr=data[pd.sigma_field], color=colors[0], label=pd.t_title) if pd.x_title and plot_bmode: @@ -607,23 +639,28 @@ def plot(self, data, colors=['r', 'b'], log_yscale=False, ax.set_xscale('log') ax.set_yscale(yscale) ax.set_xlim(xlim) + # To prevent too-long decimal y-axis ticklabels that push the label out of frame + ax.ticklabel_format(axis='y', style='sci', scilimits=(-3,5)) ax.set_ylabel(pd.y_title) ax.legend() if pd.x_field and plot_bmode and pd.t_im_field: # Both yb and y_im: plot (y, yb) on one plot and (y_im, yb_im) on the other. ax = fig.add_subplot(nrows, 1, 2) + ax.axhline(0, alpha=0.7, color='gray') ax.errorbar(data[r], data[pd.t_im_field], yerr=data[pd.sigma_field], color=colors[0], label=pd.t_im_title) ax.errorbar(data[r], data[pd.x_im_field], yerr=data[pd.sigma_field], color=colors[1], label=pd.x_im_title) ax.set_xscale('log') ax.set_yscale(yscale) + ax.ticklabel_format(axis='y', style='sci', scilimits=(-3,5)) ax.set_xlim(xlim) ax.set_ylabel(pd.y_title) ax.legend() if plot_data_only and pd.datarandom_t_field: # Plot the data-only measurements if requested curr_plot += 1 ax = fig.add_subplot(nrows, 1, 2) + ax.axhline(0, alpha=0.7, color='gray') ax.errorbar(data[r], data[pd.datarandom_t_field+'d'], yerr=data[pd.sigma_field], color=colors[0], label=pd.datarandom_t_title+'d}$') if plot_bmode and pd.datarandom_x_field: @@ -631,12 +668,14 @@ def plot(self, data, colors=['r', 'b'], log_yscale=False, color=colors[1], label=pd.datarandom_x_title+'d}$') ax.set_xscale('log') ax.set_yscale(yscale) + ax.ticklabel_format(axis='y', style='sci', scilimits=(-3,5)) ax.set_xlim(xlim) ax.set_ylabel(pd.y_title) ax.legend() # Plot the randoms-only measurements if requested if plot_random_only and pd.datarandom_t_field: ax = fig.add_subplot(nrows, 1, nrows) + ax.axhline(0, alpha=0.7, color='gray') ax.errorbar(data[r], data[pd.datarandom_t_field+'r'], yerr=data[pd.sigma_field], color=colors[0], label=pd.datarandom_t_title+'r}$') if plot_bmode and pd.datarandom_x_field: @@ -644,10 +683,12 @@ def plot(self, data, colors=['r', 'b'], log_yscale=False, color=colors[1], label=pd.datarandom_x_title+'r}$') ax.set_xscale('log') ax.set_yscale(yscale) + ax.ticklabel_format(axis='y', style='sci', scilimits=(-3,5)) ax.set_xlim(xlim) ax.set_ylabel(pd.y_title) ax.legend() ax.set_xlabel(r) + plt.tight_layout() return fig def __call__(self, *args, **kwargs): @@ -1263,7 +1304,6 @@ def __call__(self, *args, **kwargs): def getData(self): return self.data - class WhiskerPlotStarSysTest(BaseWhiskerPlotSysTest): """ A class to make WhiskerPlots of stars. @@ -2182,7 +2222,6 @@ def getStatisticsPerCCD(self, ccds, x, y, yerr=None, z=None, stat="median"): else: raise ValueError('stat should be mean or median.') - class ScatterPlotStarVsPSFG1SysTest(BaseScatterPlotSysTest): """ A class to make ScatterPlots of star vs PSF g1 values diff --git a/tests/test_correlation_functions.py b/tests/test_correlation_functions.py index cfed85b..a537768 100755 --- a/tests/test_correlation_functions.py +++ b/tests/test_correlation_functions.py @@ -45,7 +45,7 @@ def setUp(self): (0.68766, 0.68766, 0.0, 0.0, 0.0, 0.0, 0.0), (0.79877, 0.79877, 0.0, 0.0, 0.0, 0.0, 0.0), (0.92784, 0.92784, 0.0, 0.0, 0.0, 0.0, 0.0)], - dtype=[("R_nom", float), ("", float), ("", float), ("", float), + dtype=[("R_nom [deg]", float), (" [deg]", float), ("", float), ("", float), ("sigma", float), ("weight", float), ("npairs", float)]) def test_getCF(self): @@ -80,20 +80,22 @@ def test_getCF(self): realshear = stile.sys_tests.GalaxyShearSysTest() results3 = realshear(lens_data, source_data, config=stile_args) numpy.testing.assert_equal(results, results3) - + + def test_generator(self): """Make sure the CorrelationFunctionSysTest() generator returns the right objects""" - object_list = ['GalaxyShear', 'BrightStarShear', 'StarXGalaxyDensity', 'StarXGalaxyShear', + object_list = ['GalaxyShear', 'BrightStarShear', 'StarXGalaxyDensity', 'StarXGalaxyShear', 'StarXStarShear', 'GalaxyDensityCorrelation', 'StarDensityCorrelation'] for object_type in object_list: object_1 = stile.CorrelationFunctionSysTest(object_type) object_2 = eval('stile.sys_tests.'+object_type+'SysTest()') - self.assertEqual(type(object_1), type(object_2)) - - self.assertRaises(ValueError, stile.CorrelationFunctionSysTest, 'hello') - self.assertEqual(type(stile.sys_tests.BaseCorrelationFunctionSysTest()), + self.assertEqual(type(object_1),type(object_2)) + + self.assertRaises(ValueError,stile.CorrelationFunctionSysTest,'hello') + self.assertEqual(type(stile.sys_tests.BaseCorrelationFunctionSysTest()), type(stile.CorrelationFunctionSysTest())) - - -if __name__ == '__main__': + + +if __name__=='__main__': unittest.main() +