From 72ff1c2759516d20a6e1225e9c806e905cda777d Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Wed, 18 Oct 2023 06:16:09 -0700 Subject: [PATCH 1/9] Dummy commit of srv_psf_diagnostics, ignore this code --- descqa/srv_psf_diagnostics.py | 372 ++++++++++++++++++++++++++++++++++ 1 file changed, 372 insertions(+) create mode 100644 descqa/srv_psf_diagnostics.py diff --git a/descqa/srv_psf_diagnostics.py b/descqa/srv_psf_diagnostics.py new file mode 100644 index 00000000..55f6b179 --- /dev/null +++ b/descqa/srv_psf_diagnostics.py @@ -0,0 +1,372 @@ +from __future__ import print_function, division, unicode_literals, absolute_import +import os +import re +import fnmatch +from itertools import cycle +from collections import defaultdict, OrderedDict +import numpy as np +import numexpr as ne +from scipy.stats import norm, binned_statistic +import time +import sys + +from .base import BaseValidationTest, TestResult +from .plotting import plt + +if 'mpi4py' in sys.modules: + from mpi4py import MPI + from .parallel import send_to_master, get_kind + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + has_mpi = True + print('Using parallel script, invoking parallel read') +else: + size = 1 + rank = 0 + has_mpi = False + print('Using serial script') + + +__all__ = ['CheckShear','shear_from_moments'] + + + +def shear_from_moments(Ixx,Ixy,Iyy,kind='eps'): + ''' + Get shear components from second moments + ''' + if kind=='eps': + denom = Ixx + Iyy + 2.*np.sqrt(Ixx*Iyy - Ixy**2) + elif kind=='chi': + denom = Ixx + Iyy + return (Ixx-Iyy)/denom, 2*Ixy/denom + + +class PSFDiagnostics(BaseValidationTest): + """ + Make histograms of PSF values + + This makes histograms of e1 and e2 residuals, and of the fractional + PSF size excess. + """ + + def __init__(self, **kwargs): + + +class CheckEllipticity(BaseValidationTest): + """ + Check ellipticity values and second moments + """ + + def __init__(self, **kwargs): + self.catalog_filters = kwargs.get('catalog_filters', []) + self.lgndtitle_fontsize = kwargs.get('lgndtitle_fontsize', 12) + self.truncate_cat_name = kwargs.get('truncate_cat_name', False) + self.no_version = kwargs.get('no_version', False) + self.title_size = kwargs.get('title_size', 'small') + self.font_size = kwargs.get('font_size', 12) + self.legend_size = kwargs.get('legend_size', 'x-small') + self.ra = kwargs.get('ra') + self.dec = kwargs.get('dec') + self.Ixx = kwargs.get('Ixx') + self.Ixy= kwargs.get('Ixy') + self.Iyy= kwargs.get('Iyy') + self.IxxPSF= kwargs.get('IxxPSF') + self.IxyPSF= kwargs.get('IxyPSF') + self.IyyPSF= kwargs.get('IyyPSF') + self.psf_fwhm = kwargs.get('psf_fwhm') + self.bands = kwargs.get('bands') + + if not any(( + self.catalog_filters, + )): + raise ValueError('you need to specify catalog_filters for these checks, add an extendedness flag, a good flag and a magnitude range') + + self.enable_individual_summary = bool(kwargs.get('enable_individual_summary', True)) + self.enable_aggregated_summary = bool(kwargs.get('enable_aggregated_summary', False)) + self.always_show_plot = bool(kwargs.get('always_show_plot', True)) + + self.nbins = int(kwargs.get('nbins', 50)) + self.prop_cycle = None + + self.current_catalog_name = None + self.current_failed_count = None + self._aggregated_header = list() + self._aggregated_table = list() + self._individual_header = list() + self._individual_table = list() + + super(CheckEllipticity, self).__init__(**kwargs) + + def record_result(self, results, quantity_name=None, more_info=None, failed=None, individual_only=False): + if isinstance(results, dict): + self.current_failed_count += sum(1 for v in results.values() if v[1] == 'fail') + elif failed: + self.current_failed_count += 1 + + if self.enable_individual_summary: + if quantity_name is None: + self._individual_header.append(self.format_result_header(results, failed)) + else: + self._individual_table.append(self.format_result_row(results, quantity_name, more_info)) + + if self.enable_aggregated_summary and not individual_only: + if quantity_name is None: + results = '{} {}'.format(self.current_catalog_name, results) if self.current_catalog_name else results + self._aggregated_header.append(self.format_result_header(results, failed)) + else: + quantity_name = '{} {}'.format(self.current_catalog_name, quantity_name) if self.current_catalog_name else quantity_name + self._aggregated_table.append(self.format_result_row(results, quantity_name, more_info)) + + def format_result_row(self, results, quantity_name, more_info): + more_info = 'title="{}"'.format(more_info) if more_info else '' + output = ['', '{0}'.format(quantity_name, more_info)] + output.append('') + return ''.join(output) + + @staticmethod + def format_result_header(results, failed=False): + return '{0}'.format(results, 'class="fail"' if failed else '') + + def plot_moments_band(self, Ixx,Ixy,Iyy,band,output_dir): + ''' + Plot moments for each band + ''' + plt.figure() + bins = np.logspace(-1.,1.5,100) + bins_mid = (bins[1:]+bins[:-1])/2. + Ixx_out, bin_edges = np.histogram(Ixx, bins=bins) + Ixy_out, bin_edges = np.histogram(Ixy, bins=bins) + Iyy_out, bin_edges = np.histogram(Iyy, bins=bins) + self.record_result((0,'moments_'+band),'moments_'+band,'moments_'+band+'.png') + plt.plot(bins_mid,Ixx_out,'b',label='Ixx') + plt.plot(bins_mid,Iyy_out,'r--',label='Iyy') + plt.plot(bins_mid,Ixy_out,'k-.',label='Ixy') + plt.yscale('log') + plt.xscale('log') + plt.title(band) + plt.legend() + plt.savefig(os.path.join(output_dir, 'moments_'+band+'.png')) + plt.close() + return + + def plot_ellipticities_band(self,e1,e2,band,output_dir): + ''' + Plot elliptiticies for each band + ''' + plt.figure() + bins = np.linspace(0.,1.,51) + bins_mid = (bins[1:]+bins[:-1])/2. + e1_out, bin_edges = np.histogram(e1, bins=bins) + e2_out, bin_edges = np.histogram(e2, bins=bins) + self.record_result((0,'ell_'+band),'ell_'+band,'ell_'+band+'.png') + plt.plot(bins_mid,e1_out,'b',label='e1') + plt.plot(bins_mid,e2_out,'r--',label='e2') + plt.yscale('log') + #plt.xscale('log') + plt.title(band) + plt.legend() + plt.savefig(os.path.join(output_dir, 'ell_'+band+'.png')) + plt.close() + return + + def plot_psf(self,fwhm,band,output_dir): + ''' + Plot elliptiticies for each band + ''' + plt.figure() + bins = np.linspace(0.,1.5,201) + bins_mid = (bins[1:]+bins[:-1])/2. + fwhm_out, bin_edges = np.histogram(fwhm, bins=bins) + self.record_result((0,'psf_'+band),'psf_'+band,'psf_'+band+'.png') + plt.plot(bins_mid,fwhm_out,'b',label='psf') + plt.yscale('log') + plt.xlim([0.5,1.4]) + plt.title(band) + plt.legend() + plt.savefig(os.path.join(output_dir, 'psf_'+band+'.png')) + plt.close() + return + + + + + def generate_summary(self, output_dir, aggregated=False): + if aggregated: + if not self.enable_aggregated_summary: + return + header = self._aggregated_header + table = self._aggregated_table + else: + if not self.enable_individual_summary: + return + header = self._individual_header + table = self._individual_table + + with open(os.path.join(output_dir, 'SUMMARY.html'), 'w') as f: + f.write('\n') + + f.write('
\n') + + f.write('\n') + f.write('\n') + for line in table: + f.write(line) + f.write('\n') + f.write('
Quantity
\n') + + if not aggregated: + self._individual_header.clear() + self._individual_table.clear() + + def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): + + all_quantities = sorted(map(str, catalog_instance.list_all_quantities(True))) + + self.prop_cycle = cycle(iter(plt.rcParams['axes.prop_cycle'])) + self.current_catalog_name = catalog_name + self.current_failed_count = 0 + galaxy_count = None + quantity_hashes = defaultdict(set) + + if rank==0: + self.record_result('Running galaxy number test on {} {}'.format( + catalog_name, + getattr(catalog_instance, 'version', ''), + individual_only=True, + )) + + if self.truncate_cat_name: + catalog_name = catalog_name.partition("_")[0] + version = getattr(catalog_instance, 'version', '') if not self.no_version else '' + + # create filter labels + filters=[] + for i, filt in enumerate(self.catalog_filters): + filters = filt['filters'] + + print(filters) + lgnd_loc_dflt ='best' + + + label_tot=[] + plots_tot=[] + + #quantities=[self.ra,self.dec,self.Ixx,self.Iyy,self.Ixy,self.IxxPSF, self.IyyPSF, self.IxyPSF] + + # doing everything per-band first of all + for band in self.bands: + quantities=[self.Ixx+'_'+band,self.Iyy+'_'+band,self.Ixy+'_'+band,self.IxxPSF+'_'+band, self.IyyPSF+'_'+band, self.IxyPSF+'_'+band, self.psf_fwhm+'_'+band] + quantities = tuple(quantities) + + + + # reading in the data + if len(filters) > 0: + catalog_data = catalog_instance.get_quantities(quantities,filters=filters,return_iterator=False) + else: + catalog_data = catalog_instance.get_quantities(quantities,return_iterator=False) + a = time.time() + + data_rank={} + recvbuf={} + for quantity in quantities: + data_rank[quantity] = catalog_data[quantity] + if has_mpi: + if rank==0: + kind = get_kind(data_rank[quantity][0]) # assumes at least one element of data on rank 0 + else: + kind = '' + kind = comm.bcast(kind, root=0) + recvbuf[quantity] = send_to_master(data_rank[quantity],kind) + else: + recvbuf[quantity] = data_rank[quantity] + + + e1,e2 = shear_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Ixy+'_'+band],recvbuf[self.Iyy+'_'+band]) + e1psf,e2psf = shear_from_moments(recvbuf[self.IxxPSF+'_'+band],recvbuf[self.IxyPSF+'_'+band],recvbuf[self.IyyPSF+'_'+band]) + + Ixx = recvbuf[self.Ixx+'_'+band] + Iyy = recvbuf[self.Iyy+'_'+band] + Ixy = recvbuf[self.Ixy+'_'+band] + fwhm = recvbuf[self.psf_fwhm+'_'+band] + + + self.plot_moments_band(Ixx,Ixy,Iyy,band,output_dir) + self.plot_ellipticities_band(e1,e2,band,output_dir) + self.plot_psf(fwhm,band,output_dir) + + + # plot moments directly per filter. For good, star, galaxy + # FWHM of the psf + # calculate ellpiticities and make sure they're alright + # look at different bands + # note that we want to look by magnitude or SNR to understand the longer tail in moments + # PSF ellipticity whisker plot? + # look at what validate_drp is + + # s1/s2 plots + + # look at full ellipticity distribution test as well + #https://github.com/LSSTDESC/descqa/blob/master/descqa/EllipticityDistribution.py + #DC2 validation github - PSF ellipticity + # https://github.com/LSSTDESC/DC2-analysis/blob/master/validation/Run_1.2p_PSF_tests.ipynb + + #https://github.com/LSSTDESC/DC2-analysis/blob/master/validation/DC2_calexp_src_validation_1p2.ipynb + + # Look at notes here: + #https://github.com/LSSTDESC/DC2-production/issues/340 + + # get PSF FWHM directly from data, note comments on here: + # https://github.com/LSSTDESC/DC2-analysis/blob/u/wmwv/DR6_dask_refactor/validation/validate_dc2_run2.2i_object_table_dask.ipynb about focussing of the "telescope" + + + '''mask_finite = np.isfinite(e1)&np.isfinite(e2) + bs_out = bs(e1[mask_finite],values = e2[mask_finite],bins=100,statistic='mean') + plt.figure() + quantity_hashes[0].add('s1s2') + self.record_result((0,'s1s2'),'s1s2','p_s1s2.png') + plt.plot(bs_out[1][1:],bs_out[0]) + plt.savefig(os.path.join(output_dir, 'p_s1s2.png')) + plt.close() + + + plt.figure() + quantity_hashes[0].add('s1') + self.record_result((0,'s1'),'s1','p_s1.png') + #plt.hist(e1,bins=np.linspace(-1.,1.,100)) + plt.hist(e1psf,bins=100)#np.linspace(-1.,1.,100)) + plt.savefig(os.path.join(output_dir, 'p_s1.png')) + plt.close() + plt.figure() + quantity_hashes[0].add('s2') + self.record_result((0,'s2'),'s2','p_s2.png') + #plt.hist(e2,bins=np.linspace(-1.,1.,100)) + plt.hist(e2psf,bins=100)#np.linspace(-1.,1.,100)) + plt.savefig(os.path.join(output_dir, 'p_s2.png')) + plt.close()''' + '''plt.figure() + quantity_hashes[0].add('s12') + self.record_result((0,'s12'),'s12','p_s12.png') + plt.hist2d(e1,e2,bins=100) + plt.savefig(os.path.join(output_dir, 'p_s12.png')) + plt.close()''' + + if rank==0: + self.generate_summary(output_dir) + else: + self.current_failed_count=0 + + if has_mpi: + self.current_failed_count = comm.bcast(self.current_failed_count, root=0) + + return TestResult(passed=(self.current_failed_count == 0), score=self.current_failed_count) + + def conclude_test(self, output_dir): + self.generate_summary(output_dir, aggregated=True) From 0a6f078176e90a2fe87c5da0c9d6c1747847f0b9 Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Wed, 13 Dec 2023 08:25:39 -0800 Subject: [PATCH 2/9] Added new plots to srv_shear and removed srv_psf_diagnostics for now --- descqa/srv_psf_diagnostics.py | 372 ---------------------------------- descqa/srv_shear.py | 42 +++- 2 files changed, 38 insertions(+), 376 deletions(-) delete mode 100644 descqa/srv_psf_diagnostics.py diff --git a/descqa/srv_psf_diagnostics.py b/descqa/srv_psf_diagnostics.py deleted file mode 100644 index 55f6b179..00000000 --- a/descqa/srv_psf_diagnostics.py +++ /dev/null @@ -1,372 +0,0 @@ -from __future__ import print_function, division, unicode_literals, absolute_import -import os -import re -import fnmatch -from itertools import cycle -from collections import defaultdict, OrderedDict -import numpy as np -import numexpr as ne -from scipy.stats import norm, binned_statistic -import time -import sys - -from .base import BaseValidationTest, TestResult -from .plotting import plt - -if 'mpi4py' in sys.modules: - from mpi4py import MPI - from .parallel import send_to_master, get_kind - comm = MPI.COMM_WORLD - size = comm.Get_size() - rank = comm.Get_rank() - has_mpi = True - print('Using parallel script, invoking parallel read') -else: - size = 1 - rank = 0 - has_mpi = False - print('Using serial script') - - -__all__ = ['CheckShear','shear_from_moments'] - - - -def shear_from_moments(Ixx,Ixy,Iyy,kind='eps'): - ''' - Get shear components from second moments - ''' - if kind=='eps': - denom = Ixx + Iyy + 2.*np.sqrt(Ixx*Iyy - Ixy**2) - elif kind=='chi': - denom = Ixx + Iyy - return (Ixx-Iyy)/denom, 2*Ixy/denom - - -class PSFDiagnostics(BaseValidationTest): - """ - Make histograms of PSF values - - This makes histograms of e1 and e2 residuals, and of the fractional - PSF size excess. - """ - - def __init__(self, **kwargs): - - -class CheckEllipticity(BaseValidationTest): - """ - Check ellipticity values and second moments - """ - - def __init__(self, **kwargs): - self.catalog_filters = kwargs.get('catalog_filters', []) - self.lgndtitle_fontsize = kwargs.get('lgndtitle_fontsize', 12) - self.truncate_cat_name = kwargs.get('truncate_cat_name', False) - self.no_version = kwargs.get('no_version', False) - self.title_size = kwargs.get('title_size', 'small') - self.font_size = kwargs.get('font_size', 12) - self.legend_size = kwargs.get('legend_size', 'x-small') - self.ra = kwargs.get('ra') - self.dec = kwargs.get('dec') - self.Ixx = kwargs.get('Ixx') - self.Ixy= kwargs.get('Ixy') - self.Iyy= kwargs.get('Iyy') - self.IxxPSF= kwargs.get('IxxPSF') - self.IxyPSF= kwargs.get('IxyPSF') - self.IyyPSF= kwargs.get('IyyPSF') - self.psf_fwhm = kwargs.get('psf_fwhm') - self.bands = kwargs.get('bands') - - if not any(( - self.catalog_filters, - )): - raise ValueError('you need to specify catalog_filters for these checks, add an extendedness flag, a good flag and a magnitude range') - - self.enable_individual_summary = bool(kwargs.get('enable_individual_summary', True)) - self.enable_aggregated_summary = bool(kwargs.get('enable_aggregated_summary', False)) - self.always_show_plot = bool(kwargs.get('always_show_plot', True)) - - self.nbins = int(kwargs.get('nbins', 50)) - self.prop_cycle = None - - self.current_catalog_name = None - self.current_failed_count = None - self._aggregated_header = list() - self._aggregated_table = list() - self._individual_header = list() - self._individual_table = list() - - super(CheckEllipticity, self).__init__(**kwargs) - - def record_result(self, results, quantity_name=None, more_info=None, failed=None, individual_only=False): - if isinstance(results, dict): - self.current_failed_count += sum(1 for v in results.values() if v[1] == 'fail') - elif failed: - self.current_failed_count += 1 - - if self.enable_individual_summary: - if quantity_name is None: - self._individual_header.append(self.format_result_header(results, failed)) - else: - self._individual_table.append(self.format_result_row(results, quantity_name, more_info)) - - if self.enable_aggregated_summary and not individual_only: - if quantity_name is None: - results = '{} {}'.format(self.current_catalog_name, results) if self.current_catalog_name else results - self._aggregated_header.append(self.format_result_header(results, failed)) - else: - quantity_name = '{} {}'.format(self.current_catalog_name, quantity_name) if self.current_catalog_name else quantity_name - self._aggregated_table.append(self.format_result_row(results, quantity_name, more_info)) - - def format_result_row(self, results, quantity_name, more_info): - more_info = 'title="{}"'.format(more_info) if more_info else '' - output = ['', '{0}'.format(quantity_name, more_info)] - output.append('') - return ''.join(output) - - @staticmethod - def format_result_header(results, failed=False): - return '{0}'.format(results, 'class="fail"' if failed else '') - - def plot_moments_band(self, Ixx,Ixy,Iyy,band,output_dir): - ''' - Plot moments for each band - ''' - plt.figure() - bins = np.logspace(-1.,1.5,100) - bins_mid = (bins[1:]+bins[:-1])/2. - Ixx_out, bin_edges = np.histogram(Ixx, bins=bins) - Ixy_out, bin_edges = np.histogram(Ixy, bins=bins) - Iyy_out, bin_edges = np.histogram(Iyy, bins=bins) - self.record_result((0,'moments_'+band),'moments_'+band,'moments_'+band+'.png') - plt.plot(bins_mid,Ixx_out,'b',label='Ixx') - plt.plot(bins_mid,Iyy_out,'r--',label='Iyy') - plt.plot(bins_mid,Ixy_out,'k-.',label='Ixy') - plt.yscale('log') - plt.xscale('log') - plt.title(band) - plt.legend() - plt.savefig(os.path.join(output_dir, 'moments_'+band+'.png')) - plt.close() - return - - def plot_ellipticities_band(self,e1,e2,band,output_dir): - ''' - Plot elliptiticies for each band - ''' - plt.figure() - bins = np.linspace(0.,1.,51) - bins_mid = (bins[1:]+bins[:-1])/2. - e1_out, bin_edges = np.histogram(e1, bins=bins) - e2_out, bin_edges = np.histogram(e2, bins=bins) - self.record_result((0,'ell_'+band),'ell_'+band,'ell_'+band+'.png') - plt.plot(bins_mid,e1_out,'b',label='e1') - plt.plot(bins_mid,e2_out,'r--',label='e2') - plt.yscale('log') - #plt.xscale('log') - plt.title(band) - plt.legend() - plt.savefig(os.path.join(output_dir, 'ell_'+band+'.png')) - plt.close() - return - - def plot_psf(self,fwhm,band,output_dir): - ''' - Plot elliptiticies for each band - ''' - plt.figure() - bins = np.linspace(0.,1.5,201) - bins_mid = (bins[1:]+bins[:-1])/2. - fwhm_out, bin_edges = np.histogram(fwhm, bins=bins) - self.record_result((0,'psf_'+band),'psf_'+band,'psf_'+band+'.png') - plt.plot(bins_mid,fwhm_out,'b',label='psf') - plt.yscale('log') - plt.xlim([0.5,1.4]) - plt.title(band) - plt.legend() - plt.savefig(os.path.join(output_dir, 'psf_'+band+'.png')) - plt.close() - return - - - - - def generate_summary(self, output_dir, aggregated=False): - if aggregated: - if not self.enable_aggregated_summary: - return - header = self._aggregated_header - table = self._aggregated_table - else: - if not self.enable_individual_summary: - return - header = self._individual_header - table = self._individual_table - - with open(os.path.join(output_dir, 'SUMMARY.html'), 'w') as f: - f.write('\n') - - f.write('
\n') - - f.write('\n') - f.write('\n') - for line in table: - f.write(line) - f.write('\n') - f.write('
Quantity
\n') - - if not aggregated: - self._individual_header.clear() - self._individual_table.clear() - - def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): - - all_quantities = sorted(map(str, catalog_instance.list_all_quantities(True))) - - self.prop_cycle = cycle(iter(plt.rcParams['axes.prop_cycle'])) - self.current_catalog_name = catalog_name - self.current_failed_count = 0 - galaxy_count = None - quantity_hashes = defaultdict(set) - - if rank==0: - self.record_result('Running galaxy number test on {} {}'.format( - catalog_name, - getattr(catalog_instance, 'version', ''), - individual_only=True, - )) - - if self.truncate_cat_name: - catalog_name = catalog_name.partition("_")[0] - version = getattr(catalog_instance, 'version', '') if not self.no_version else '' - - # create filter labels - filters=[] - for i, filt in enumerate(self.catalog_filters): - filters = filt['filters'] - - print(filters) - lgnd_loc_dflt ='best' - - - label_tot=[] - plots_tot=[] - - #quantities=[self.ra,self.dec,self.Ixx,self.Iyy,self.Ixy,self.IxxPSF, self.IyyPSF, self.IxyPSF] - - # doing everything per-band first of all - for band in self.bands: - quantities=[self.Ixx+'_'+band,self.Iyy+'_'+band,self.Ixy+'_'+band,self.IxxPSF+'_'+band, self.IyyPSF+'_'+band, self.IxyPSF+'_'+band, self.psf_fwhm+'_'+band] - quantities = tuple(quantities) - - - - # reading in the data - if len(filters) > 0: - catalog_data = catalog_instance.get_quantities(quantities,filters=filters,return_iterator=False) - else: - catalog_data = catalog_instance.get_quantities(quantities,return_iterator=False) - a = time.time() - - data_rank={} - recvbuf={} - for quantity in quantities: - data_rank[quantity] = catalog_data[quantity] - if has_mpi: - if rank==0: - kind = get_kind(data_rank[quantity][0]) # assumes at least one element of data on rank 0 - else: - kind = '' - kind = comm.bcast(kind, root=0) - recvbuf[quantity] = send_to_master(data_rank[quantity],kind) - else: - recvbuf[quantity] = data_rank[quantity] - - - e1,e2 = shear_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Ixy+'_'+band],recvbuf[self.Iyy+'_'+band]) - e1psf,e2psf = shear_from_moments(recvbuf[self.IxxPSF+'_'+band],recvbuf[self.IxyPSF+'_'+band],recvbuf[self.IyyPSF+'_'+band]) - - Ixx = recvbuf[self.Ixx+'_'+band] - Iyy = recvbuf[self.Iyy+'_'+band] - Ixy = recvbuf[self.Ixy+'_'+band] - fwhm = recvbuf[self.psf_fwhm+'_'+band] - - - self.plot_moments_band(Ixx,Ixy,Iyy,band,output_dir) - self.plot_ellipticities_band(e1,e2,band,output_dir) - self.plot_psf(fwhm,band,output_dir) - - - # plot moments directly per filter. For good, star, galaxy - # FWHM of the psf - # calculate ellpiticities and make sure they're alright - # look at different bands - # note that we want to look by magnitude or SNR to understand the longer tail in moments - # PSF ellipticity whisker plot? - # look at what validate_drp is - - # s1/s2 plots - - # look at full ellipticity distribution test as well - #https://github.com/LSSTDESC/descqa/blob/master/descqa/EllipticityDistribution.py - #DC2 validation github - PSF ellipticity - # https://github.com/LSSTDESC/DC2-analysis/blob/master/validation/Run_1.2p_PSF_tests.ipynb - - #https://github.com/LSSTDESC/DC2-analysis/blob/master/validation/DC2_calexp_src_validation_1p2.ipynb - - # Look at notes here: - #https://github.com/LSSTDESC/DC2-production/issues/340 - - # get PSF FWHM directly from data, note comments on here: - # https://github.com/LSSTDESC/DC2-analysis/blob/u/wmwv/DR6_dask_refactor/validation/validate_dc2_run2.2i_object_table_dask.ipynb about focussing of the "telescope" - - - '''mask_finite = np.isfinite(e1)&np.isfinite(e2) - bs_out = bs(e1[mask_finite],values = e2[mask_finite],bins=100,statistic='mean') - plt.figure() - quantity_hashes[0].add('s1s2') - self.record_result((0,'s1s2'),'s1s2','p_s1s2.png') - plt.plot(bs_out[1][1:],bs_out[0]) - plt.savefig(os.path.join(output_dir, 'p_s1s2.png')) - plt.close() - - - plt.figure() - quantity_hashes[0].add('s1') - self.record_result((0,'s1'),'s1','p_s1.png') - #plt.hist(e1,bins=np.linspace(-1.,1.,100)) - plt.hist(e1psf,bins=100)#np.linspace(-1.,1.,100)) - plt.savefig(os.path.join(output_dir, 'p_s1.png')) - plt.close() - plt.figure() - quantity_hashes[0].add('s2') - self.record_result((0,'s2'),'s2','p_s2.png') - #plt.hist(e2,bins=np.linspace(-1.,1.,100)) - plt.hist(e2psf,bins=100)#np.linspace(-1.,1.,100)) - plt.savefig(os.path.join(output_dir, 'p_s2.png')) - plt.close()''' - '''plt.figure() - quantity_hashes[0].add('s12') - self.record_result((0,'s12'),'s12','p_s12.png') - plt.hist2d(e1,e2,bins=100) - plt.savefig(os.path.join(output_dir, 'p_s12.png')) - plt.close()''' - - if rank==0: - self.generate_summary(output_dir) - else: - self.current_failed_count=0 - - if has_mpi: - self.current_failed_count = comm.bcast(self.current_failed_count, root=0) - - return TestResult(passed=(self.current_failed_count == 0), score=self.current_failed_count) - - def conclude_test(self, output_dir): - self.generate_summary(output_dir, aggregated=True) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index 8fec7f8c..b0fa3f5c 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -42,6 +42,8 @@ def shear_from_moments(Ixx,Ixy,Iyy,kind='eps'): denom = Ixx + Iyy return (Ixx-Iyy)/denom, 2*Ixy/denom +def size_from_moments(Ixx, Iyy): + return Ixx + Iyy class CheckEllipticity(BaseValidationTest): """ @@ -162,7 +164,7 @@ def plot_ellipticities_band(self,e1,e2,band,output_dir): def plot_psf(self,fwhm,band,output_dir): ''' - Plot elliptiticies for each band + Plot PSF for each band ''' plt.figure() bins = np.linspace(0.,1.5,201) @@ -178,8 +180,37 @@ def plot_psf(self,fwhm,band,output_dir): plt.close() return - - + def plot_e1e2_residuals(self,e1,e2,epsf1,epsf2): + ''' + Plot e1,e2 residuals with respect to model + ''' + plt.figure() + bins = np.linspace(-1.,1.,201) + bins_mid = (bins[1:]+bins[:-1])/2. + e1residual_out, bin_edges = np.histogram(e1-epsf1, bins=bins) + e2residual_out, bin_edges = np.histogram(e2-epsf2, bins=bins) + self.record_result((0,'e1e2_residuals_'+band),'e1e2_residuals_'+band,'e1e2_residuals_'+band+'.png') + plt.plot(bins_mid,e1residual_out,'b',label='e1-e1psf') + plt.plot(bins_mid,e2residual_out,'r--',label='e2-e2psf') + plt.title('e1/e2 residuals vs model, band'+band) + plt.savefig(os.path.join(output_dir, 'e1e2_residuals_'+band+'.png')) + plt.close() + return + + def plot_Tfrac_residuals(self,T,Tpsf): + ''' + Plot T fractional residuals with respect to model + ''' + plt.figure() + bins = np.linspace(-0.1,0.1,201) + bins_mid = (bins[1:]+bins[:-1])/2. + tresidual_out, bin_edges = np.histogram((T-Tpsf)/Tpsf, bins=bins) + self.record_result((0,'tfrac_residuals_'+band),'tfrac_residuals_'+band,'tfrac_residuals_'+band+'.png') + plt.plot(bins_mid,tresidual_out,'b',label='(T-Tpsf)/Tpsf') + plt.title('T fractional residuals vs model, band'+band) + plt.savefig(os.path.join(output_dir, 'tfrac_residuals_'+band+'.png')) + plt.close() + return def generate_summary(self, output_dir, aggregated=False): if aggregated: @@ -280,16 +311,19 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): e1,e2 = shear_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Ixy+'_'+band],recvbuf[self.Iyy+'_'+band]) e1psf,e2psf = shear_from_moments(recvbuf[self.IxxPSF+'_'+band],recvbuf[self.IxyPSF+'_'+band],recvbuf[self.IyyPSF+'_'+band]) + T = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) + Tpsf = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) Ixx = recvbuf[self.Ixx+'_'+band] Iyy = recvbuf[self.Iyy+'_'+band] Ixy = recvbuf[self.Ixy+'_'+band] fwhm = recvbuf[self.psf_fwhm+'_'+band] - self.plot_moments_band(Ixx,Ixy,Iyy,band,output_dir) self.plot_ellipticities_band(e1,e2,band,output_dir) self.plot_psf(fwhm,band,output_dir) + self.plot_e_residuals(e1,e2,e1psf,e2psf) + self.plot_Tfrac_residuals(T,Tpsf) # plot moments directly per filter. For good, star, galaxy From 4e5a9bd010d15f9c0a0f8de667d26f88fd8acc58 Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Thu, 14 Dec 2023 06:37:00 -0800 Subject: [PATCH 3/9] Fixed variable names in srv_shear --- descqa/srv_shear.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index b0fa3f5c..bd40476a 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -180,7 +180,7 @@ def plot_psf(self,fwhm,band,output_dir): plt.close() return - def plot_e1e2_residuals(self,e1,e2,epsf1,epsf2): + def plot_e1e2_residuals(self,e1,e2,epsf1,epsf2,band,output_dir): ''' Plot e1,e2 residuals with respect to model ''' @@ -197,7 +197,7 @@ def plot_e1e2_residuals(self,e1,e2,epsf1,epsf2): plt.close() return - def plot_Tfrac_residuals(self,T,Tpsf): + def plot_Tfrac_residuals(self,T,Tpsf,band,output_dir): ''' Plot T fractional residuals with respect to model ''' @@ -322,8 +322,8 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): self.plot_moments_band(Ixx,Ixy,Iyy,band,output_dir) self.plot_ellipticities_band(e1,e2,band,output_dir) self.plot_psf(fwhm,band,output_dir) - self.plot_e_residuals(e1,e2,e1psf,e2psf) - self.plot_Tfrac_residuals(T,Tpsf) + self.plot_e1e2_residuals(e1,e2,e1psf,e2psf,band,output_dir) + self.plot_Tfrac_residuals(T,Tpsf,band,output_dir) # plot moments directly per filter. For good, star, galaxy From 7f89655d053edbff461442a18442ec29501ff526 Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Thu, 14 Dec 2023 06:49:41 -0800 Subject: [PATCH 4/9] Eliminating a lot of old comments in srv_shear --- descqa/srv_shear.py | 62 ++------------------------------------------- 1 file changed, 2 insertions(+), 60 deletions(-) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index bd40476a..cc698a0e 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -278,15 +278,11 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): label_tot=[] plots_tot=[] - #quantities=[self.ra,self.dec,self.Ixx,self.Iyy,self.Ixy,self.IxxPSF, self.IyyPSF, self.IxyPSF] - # doing everything per-band first of all for band in self.bands: quantities=[self.Ixx+'_'+band,self.Iyy+'_'+band,self.Ixy+'_'+band,self.IxxPSF+'_'+band, self.IyyPSF+'_'+band, self.IxyPSF+'_'+band, self.psf_fwhm+'_'+band] quantities = tuple(quantities) - - # reading in the data if len(filters) > 0: catalog_data = catalog_instance.get_quantities(quantities,filters=filters,return_iterator=False) @@ -313,6 +309,8 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): e1psf,e2psf = shear_from_moments(recvbuf[self.IxxPSF+'_'+band],recvbuf[self.IxyPSF+'_'+band],recvbuf[self.IyyPSF+'_'+band]) T = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) Tpsf = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) + de1 = e1-e1psf + de2 = e2-e2psf Ixx = recvbuf[self.Ixx+'_'+band] Iyy = recvbuf[self.Iyy+'_'+band] @@ -325,62 +323,6 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): self.plot_e1e2_residuals(e1,e2,e1psf,e2psf,band,output_dir) self.plot_Tfrac_residuals(T,Tpsf,band,output_dir) - - # plot moments directly per filter. For good, star, galaxy - # FWHM of the psf - # calculate ellpiticities and make sure they're alright - # look at different bands - # note that we want to look by magnitude or SNR to understand the longer tail in moments - # PSF ellipticity whisker plot? - # look at what validate_drp is - - # s1/s2 plots - - # look at full ellipticity distribution test as well - #https://github.com/LSSTDESC/descqa/blob/master/descqa/EllipticityDistribution.py - #DC2 validation github - PSF ellipticity - # https://github.com/LSSTDESC/DC2-analysis/blob/master/validation/Run_1.2p_PSF_tests.ipynb - - #https://github.com/LSSTDESC/DC2-analysis/blob/master/validation/DC2_calexp_src_validation_1p2.ipynb - - # Look at notes here: - #https://github.com/LSSTDESC/DC2-production/issues/340 - - # get PSF FWHM directly from data, note comments on here: - # https://github.com/LSSTDESC/DC2-analysis/blob/u/wmwv/DR6_dask_refactor/validation/validate_dc2_run2.2i_object_table_dask.ipynb about focussing of the "telescope" - - - '''mask_finite = np.isfinite(e1)&np.isfinite(e2) - bs_out = bs(e1[mask_finite],values = e2[mask_finite],bins=100,statistic='mean') - plt.figure() - quantity_hashes[0].add('s1s2') - self.record_result((0,'s1s2'),'s1s2','p_s1s2.png') - plt.plot(bs_out[1][1:],bs_out[0]) - plt.savefig(os.path.join(output_dir, 'p_s1s2.png')) - plt.close() - - - plt.figure() - quantity_hashes[0].add('s1') - self.record_result((0,'s1'),'s1','p_s1.png') - #plt.hist(e1,bins=np.linspace(-1.,1.,100)) - plt.hist(e1psf,bins=100)#np.linspace(-1.,1.,100)) - plt.savefig(os.path.join(output_dir, 'p_s1.png')) - plt.close() - plt.figure() - quantity_hashes[0].add('s2') - self.record_result((0,'s2'),'s2','p_s2.png') - #plt.hist(e2,bins=np.linspace(-1.,1.,100)) - plt.hist(e2psf,bins=100)#np.linspace(-1.,1.,100)) - plt.savefig(os.path.join(output_dir, 'p_s2.png')) - plt.close()''' - '''plt.figure() - quantity_hashes[0].add('s12') - self.record_result((0,'s12'),'s12','p_s12.png') - plt.hist2d(e1,e2,bins=100) - plt.savefig(os.path.join(output_dir, 'p_s12.png')) - plt.close()''' - if rank==0: self.generate_summary(output_dir) else: From 1c45dbf49f6828636e3b1ad6c4a533176dc3293a Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Thu, 4 Jan 2024 06:10:05 -0800 Subject: [PATCH 5/9] Minor indentation fix --- descqa/srv_shear.py | 87 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index cc698a0e..bfe7ef90 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -9,9 +9,12 @@ from scipy.stats import norm, binned_statistic import time import sys +import treecorr from .base import BaseValidationTest, TestResult from .plotting import plt +import matplotlib.transforms as mtrans + if 'mpi4py' in sys.modules: from mpi4py import MPI @@ -45,6 +48,20 @@ def shear_from_moments(Ixx,Ixy,Iyy,kind='eps'): def size_from_moments(Ixx, Iyy): return Ixx + Iyy +def compute_rowe(self, i, ra, dec, q1, q2): + n = len(ra) + print(f"Computing Rowe statistic rho_{i} from {n} objects") + + corr = treecorr.GGCorrelation(self.config) + cat1 = treecorr.Catalog( + ra=ra, dec=dec, g1=q1[0], g2=q1[1], ra_units="deg", dec_units="deg" + ) + cat2 = treecorr.Catalog( + ra=ra, dec=dec, g1=q2[0], g2=q2[1], ra_units="deg", dec_units="deg" + ) + corr.process(cat1, cat2) + return corr.meanr, corr.xip, corr.varxip**0.5 + class CheckEllipticity(BaseValidationTest): """ Check ellipticity values and second moments @@ -68,6 +85,7 @@ def __init__(self, **kwargs): self.IyyPSF= kwargs.get('IyyPSF') self.psf_fwhm = kwargs.get('psf_fwhm') self.bands = kwargs.get('bands') + self.compute_rowe = kwargs.get('rowe', False) if not any(( self.catalog_filters, @@ -211,6 +229,60 @@ def plot_Tfrac_residuals(self,T,Tpsf,band,output_dir): plt.savefig(os.path.join(output_dir, 'tfrac_residuals_'+band+'.png')) plt.close() return + + def plot_rowe_stats(self,rowe_stats,output_dir): + ''' + Plot Rowe coefficients + ''' + filename = output_dir + "rowe134.png" + ax = plt.subplot(1, 1, 1) + for j, i in enumerate([1, 3, 4]): + theta, xi, err = rowe_stats[i] + tr = mtrans.offset_copy( + ax.transData, filename, 0.05 * (j - 1), 0, units="inches" + ) + plt.errorbar( + theta, + abs(xi), + err, + fmt=".", + label=rf"$\rho_{i}$", + capsize=3, + transform=tr, + ) + plt.bar(0.0, 2e-05, width=5, align="edge", color="gray", alpha=0.2) + plt.bar(5, 1e-07, width=245, align="edge", color="gray", alpha=0.2) + plt.xscale("log") + plt.yscale("log") + plt.xlabel(r"$\theta$") + plt.ylabel(r"$\xi_+(\theta)$") + plt.legend() + + filename = output_dir + "rowe25.png" + ax = plt.subplot(1, 1, 1) + for j, i in enumerate([2, 5]): + theta, xi, err = rowe_stats[i] + tr = mtrans.offset_copy( + ax.transData, filename, 0.05 * j - 0.025, 0, units="inches" + ) + plt.errorbar( + theta, + abs(xi), + err, + fmt=".", + label=rf"$\rho_{i}$", + capsize=3, + transform=tr, + ) + plt.bar(0.0, 2e-05, width=5, align="edge", color="gray", alpha=0.2) + plt.bar(5, 1e-07, width=245, align="edge", color="gray", alpha=0.2) + plt.xscale("log") + plt.yscale("log") + plt.xlabel(r"$\theta$") + plt.ylabel(r"$\xi_+(\theta)$") + plt.legend() + + return def generate_summary(self, output_dir, aggregated=False): if aggregated: @@ -309,6 +381,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): e1psf,e2psf = shear_from_moments(recvbuf[self.IxxPSF+'_'+band],recvbuf[self.IxyPSF+'_'+band],recvbuf[self.IyyPSF+'_'+band]) T = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) Tpsf = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) + T_f = (T-Tpsf)/Tpsf de1 = e1-e1psf de2 = e2-e2psf @@ -322,6 +395,20 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): self.plot_psf(fwhm,band,output_dir) self.plot_e1e2_residuals(e1,e2,e1psf,e2psf,band,output_dir) self.plot_Tfrac_residuals(T,Tpsf,band,output_dir) + + rowe_stats = np.empty([6]) + rowe_stats[0] = 0. #dummy value, so that in the rest of the array the index + #coincides with the usual Rowe number + + if self.compute_rowe: + rowe_stats[1] = self.compute_rowe(1, ra, dec, (de1,de2), (de1,de2)) + rowe_stats[2] = self.compute_rowe(2, ra, dec, (e1,e2), (de1,de2)) + rowe_stats[3] = self.compute_rowe(3, ra, dec, (e1,e2) * T_f, (e1,e2) * T_f) + rowe_stats[4] = self.compute_rowe(4, ra, dec, (de1,de2), (e1,e2) * T_f) + rowe_stats[5] = self.compute_rowe(5, ra, dec, (e1,e2), (e1,e2) * T_f) + + self.plot_rowe_stats(rowe_stats,output_dir) + if rank==0: self.generate_summary(output_dir) From c3f7404c10ffddd09df5fd211bab8b2d9a5e57b1 Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Thu, 4 Jan 2024 06:18:43 -0800 Subject: [PATCH 6/9] Dummy changes, ignore --- descqa/srv_shear.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index bfe7ef90..ec1804ee 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -398,7 +398,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): rowe_stats = np.empty([6]) rowe_stats[0] = 0. #dummy value, so that in the rest of the array the index - #coincides with the usual Rowe number + #coincides with the usual Rowe number if self.compute_rowe: rowe_stats[1] = self.compute_rowe(1, ra, dec, (de1,de2), (de1,de2)) From a6dea9886583215674cda1aea0cfa354edae477e Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Mon, 18 Mar 2024 15:31:20 -0700 Subject: [PATCH 7/9] added Rowe statistics --- descqa/srv_shear.py | 89 ++++++++++++++++++++++++++++++--------------- 1 file changed, 59 insertions(+), 30 deletions(-) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index ec1804ee..6f1139bb 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -48,19 +48,6 @@ def shear_from_moments(Ixx,Ixy,Iyy,kind='eps'): def size_from_moments(Ixx, Iyy): return Ixx + Iyy -def compute_rowe(self, i, ra, dec, q1, q2): - n = len(ra) - print(f"Computing Rowe statistic rho_{i} from {n} objects") - - corr = treecorr.GGCorrelation(self.config) - cat1 = treecorr.Catalog( - ra=ra, dec=dec, g1=q1[0], g2=q1[1], ra_units="deg", dec_units="deg" - ) - cat2 = treecorr.Catalog( - ra=ra, dec=dec, g1=q2[0], g2=q2[1], ra_units="deg", dec_units="deg" - ) - corr.process(cat1, cat2) - return corr.meanr, corr.xip, corr.varxip**0.5 class CheckEllipticity(BaseValidationTest): """ @@ -85,7 +72,16 @@ def __init__(self, **kwargs): self.IyyPSF= kwargs.get('IyyPSF') self.psf_fwhm = kwargs.get('psf_fwhm') self.bands = kwargs.get('bands') - self.compute_rowe = kwargs.get('rowe', False) + self.compute_rowe_flag = kwargs.get('rowe', True) + self.treecorr_config = { + "min_sep": 0.5, + "max_sep": 250.0, + "nbins": 20, + "bin_slop": 0.01, + "sep_units": "arcmin", + "psf_size_units": "sigma", + } + if not any(( self.catalog_filters, @@ -236,10 +232,11 @@ def plot_rowe_stats(self,rowe_stats,output_dir): ''' filename = output_dir + "rowe134.png" ax = plt.subplot(1, 1, 1) + fig = plt.figure() for j, i in enumerate([1, 3, 4]): theta, xi, err = rowe_stats[i] tr = mtrans.offset_copy( - ax.transData, filename, 0.05 * (j - 1), 0, units="inches" + ax.transData, fig, 0.05 * (j - 1), 0, units="inches" ) plt.errorbar( theta, @@ -248,8 +245,9 @@ def plot_rowe_stats(self,rowe_stats,output_dir): fmt=".", label=rf"$\rho_{i}$", capsize=3, - transform=tr, ) + # transform=tr, + #) plt.bar(0.0, 2e-05, width=5, align="edge", color="gray", alpha=0.2) plt.bar(5, 1e-07, width=245, align="edge", color="gray", alpha=0.2) plt.xscale("log") @@ -257,13 +255,16 @@ def plot_rowe_stats(self,rowe_stats,output_dir): plt.xlabel(r"$\theta$") plt.ylabel(r"$\xi_+(\theta)$") plt.legend() + plt.savefig(filename) + plt.close() filename = output_dir + "rowe25.png" ax = plt.subplot(1, 1, 1) + fig = plt.figure() for j, i in enumerate([2, 5]): theta, xi, err = rowe_stats[i] tr = mtrans.offset_copy( - ax.transData, filename, 0.05 * j - 0.025, 0, units="inches" + ax.transData, fig, 0.05 * j - 0.025, 0, units="inches" ) plt.errorbar( theta, @@ -272,8 +273,9 @@ def plot_rowe_stats(self,rowe_stats,output_dir): fmt=".", label=rf"$\rho_{i}$", capsize=3, - transform=tr, ) + # transform=tr, + #) plt.bar(0.0, 2e-05, width=5, align="edge", color="gray", alpha=0.2) plt.bar(5, 1e-07, width=245, align="edge", color="gray", alpha=0.2) plt.xscale("log") @@ -281,7 +283,9 @@ def plot_rowe_stats(self,rowe_stats,output_dir): plt.xlabel(r"$\theta$") plt.ylabel(r"$\xi_+(\theta)$") plt.legend() - + plt.savefig(filename) + plt.close() + return def generate_summary(self, output_dir, aggregated=False): @@ -317,6 +321,16 @@ def generate_summary(self, output_dir, aggregated=False): self._individual_header.clear() self._individual_table.clear() + def compute_rowe(self, i, ra, dec, q1, q2): + n = len(ra) + print(f"Computing Rowe statistic rho_{i} from {n} objects") + + corr = treecorr.GGCorrelation(self.treecorr_config) + cat1 = treecorr.Catalog(ra=np.array(ra), dec=np.array(dec), g1=np.nan_to_num(q1[0], copy=True, nan=0.0, posinf=None, neginf=None), g2=np.nan_to_num(q1[1], copy=True, nan=0.0, posinf=None, neginf=None), ra_units="deg", dec_units="deg") + cat2 = treecorr.Catalog(ra=np.array(ra), dec=np.array(dec), g1=np.nan_to_num(q2[0], copy=True, nan=0.0, posinf=None, neginf=None), g2=np.nan_to_num(q2[1], copy=True, nan=0.0, posinf=None, neginf=None), ra_units="deg", dec_units="deg") + corr.process(cat1, cat2) #error is here + return corr.meanr, corr.xip, corr.varxip**0.5 + def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): all_quantities = sorted(map(str, catalog_instance.list_all_quantities(True))) @@ -352,7 +366,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): # doing everything per-band first of all for band in self.bands: - quantities=[self.Ixx+'_'+band,self.Iyy+'_'+band,self.Ixy+'_'+band,self.IxxPSF+'_'+band, self.IyyPSF+'_'+band, self.IxyPSF+'_'+band, self.psf_fwhm+'_'+band] + quantities=[self.Ixx+'_'+band,self.Iyy+'_'+band,self.Ixy+'_'+band,self.IxxPSF+'_'+band, self.IyyPSF+'_'+band, self.IxyPSF+'_'+band, self.psf_fwhm+'_'+band, self.ra, self.dec] quantities = tuple(quantities) # reading in the data @@ -380,7 +394,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): e1,e2 = shear_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Ixy+'_'+band],recvbuf[self.Iyy+'_'+band]) e1psf,e2psf = shear_from_moments(recvbuf[self.IxxPSF+'_'+band],recvbuf[self.IxyPSF+'_'+band],recvbuf[self.IyyPSF+'_'+band]) T = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) - Tpsf = size_from_moments(recvbuf[self.Ixx+'_'+band],recvbuf[self.Iyy+'_'+band]) + Tpsf = size_from_moments(recvbuf[self.IxxPSF+'_'+band],recvbuf[self.IyyPSF+'_'+band]) T_f = (T-Tpsf)/Tpsf de1 = e1-e1psf de2 = e2-e2psf @@ -396,18 +410,33 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): self.plot_e1e2_residuals(e1,e2,e1psf,e2psf,band,output_dir) self.plot_Tfrac_residuals(T,Tpsf,band,output_dir) - rowe_stats = np.empty([6]) - rowe_stats[0] = 0. #dummy value, so that in the rest of the array the index + #rowe_stats = np.empty([6]) + #rowe_stats[0] = 0. #dummy value, so that in the rest of the array the index #coincides with the usual Rowe number + rowe_stats = {} + + print(type(recvbuf)) + print(recvbuf[self.ra]) - if self.compute_rowe: - rowe_stats[1] = self.compute_rowe(1, ra, dec, (de1,de2), (de1,de2)) - rowe_stats[2] = self.compute_rowe(2, ra, dec, (e1,e2), (de1,de2)) - rowe_stats[3] = self.compute_rowe(3, ra, dec, (e1,e2) * T_f, (e1,e2) * T_f) - rowe_stats[4] = self.compute_rowe(4, ra, dec, (de1,de2), (e1,e2) * T_f) - rowe_stats[5] = self.compute_rowe(5, ra, dec, (e1,e2), (e1,e2) * T_f) + if self.compute_rowe_flag: + print("Computing Rowe_1") + rowe_stats[1] = self.compute_rowe(1, recvbuf[self.ra], recvbuf[self.dec], (de1,de2), (de1,de2)) + print(rowe_stats[1]) + print("Computing Rowe_2") + rowe_stats[2] = self.compute_rowe(2, recvbuf[self.ra], recvbuf[self.dec], (e1,e2), (de1,de2)) + print(rowe_stats[2]) + print("Computing Rowe_3") + rowe_stats[3] = self.compute_rowe(3, recvbuf[self.ra], recvbuf[self.dec], (e1,e2) * T_f, (e1,e2) * T_f) + print(rowe_stats[3]) + print("Computing Rowe_4") + rowe_stats[4] = self.compute_rowe(4, recvbuf[self.ra], recvbuf[self.dec], (de1,de2), (e1,e2) * T_f) + print(rowe_stats[4]) + print("Computing Rowe_5") + rowe_stats[5] = self.compute_rowe(5, recvbuf[self.ra], recvbuf[self.dec], (e1,e2), (e1,e2) * T_f) + print(rowe_stats[5]) - self.plot_rowe_stats(rowe_stats,output_dir) + print("Now plotting") + self.plot_rowe_stats(rowe_stats,output_dir) if rank==0: From 744c7967da971a85c8f59d7fd576852d89267492 Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Fri, 26 Apr 2024 05:48:55 -0700 Subject: [PATCH 8/9] Adding legends to data vs model plots --- descqa/srv_shear.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index 6f1139bb..8ea6a46d 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -206,8 +206,9 @@ def plot_e1e2_residuals(self,e1,e2,epsf1,epsf2,band,output_dir): self.record_result((0,'e1e2_residuals_'+band),'e1e2_residuals_'+band,'e1e2_residuals_'+band+'.png') plt.plot(bins_mid,e1residual_out,'b',label='e1-e1psf') plt.plot(bins_mid,e2residual_out,'r--',label='e2-e2psf') - plt.title('e1/e2 residuals vs model, band'+band) + plt.title('e1/e2 residuals vs model, band '+band) plt.savefig(os.path.join(output_dir, 'e1e2_residuals_'+band+'.png')) + plt.legend() plt.close() return @@ -221,8 +222,9 @@ def plot_Tfrac_residuals(self,T,Tpsf,band,output_dir): tresidual_out, bin_edges = np.histogram((T-Tpsf)/Tpsf, bins=bins) self.record_result((0,'tfrac_residuals_'+band),'tfrac_residuals_'+band,'tfrac_residuals_'+band+'.png') plt.plot(bins_mid,tresidual_out,'b',label='(T-Tpsf)/Tpsf') - plt.title('T fractional residuals vs model, band'+band) + plt.title('T fractional residuals vs model, band '+band) plt.savefig(os.path.join(output_dir, 'tfrac_residuals_'+band+'.png')) + plt.legend() plt.close() return From 09028142e38e00e7851d61ac31779a0598b60f41 Mon Sep 17 00:00:00 2001 From: Nacho Sevilla Date: Fri, 26 Apr 2024 06:39:06 -0700 Subject: [PATCH 9/9] Bug fix on legend plotting --- descqa/srv_shear.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/descqa/srv_shear.py b/descqa/srv_shear.py index 8ea6a46d..eef7f8d3 100644 --- a/descqa/srv_shear.py +++ b/descqa/srv_shear.py @@ -207,8 +207,8 @@ def plot_e1e2_residuals(self,e1,e2,epsf1,epsf2,band,output_dir): plt.plot(bins_mid,e1residual_out,'b',label='e1-e1psf') plt.plot(bins_mid,e2residual_out,'r--',label='e2-e2psf') plt.title('e1/e2 residuals vs model, band '+band) - plt.savefig(os.path.join(output_dir, 'e1e2_residuals_'+band+'.png')) plt.legend() + plt.savefig(os.path.join(output_dir, 'e1e2_residuals_'+band+'.png')) plt.close() return @@ -223,8 +223,8 @@ def plot_Tfrac_residuals(self,T,Tpsf,band,output_dir): self.record_result((0,'tfrac_residuals_'+band),'tfrac_residuals_'+band,'tfrac_residuals_'+band+'.png') plt.plot(bins_mid,tresidual_out,'b',label='(T-Tpsf)/Tpsf') plt.title('T fractional residuals vs model, band '+band) - plt.savefig(os.path.join(output_dir, 'tfrac_residuals_'+band+'.png')) plt.legend() + plt.savefig(os.path.join(output_dir, 'tfrac_residuals_'+band+'.png')) plt.close() return