diff --git a/descqa/AGN_NumberDensities.py b/descqa/AGN_NumberDensities.py new file mode 100644 index 00000000..b2f5b031 --- /dev/null +++ b/descqa/AGN_NumberDensities.py @@ -0,0 +1,414 @@ +from __future__ import unicode_literals, absolute_import, division +import os +import re +import numpy as np +from scipy.interpolate import interp1d +from GCR import GCRQuery +try: + from itertools import zip_longest +except ImportError: + from itertools import izip_longest as zip_longest +from .plotting import plt +from .base import BaseValidationTest, TestResult + +possible_observations = { + 'SDSS': { + 'filename_template': 'AGNdata/SDSS/richards_2006_table2.dat', + 'usecols': {'g': {'N {}'.format(self.z_lim[0]), 'redshift < {}'.format(self.z_lim[1])] + Mag_filters = ['{} < {}'.format(Mag_field_key, self.Mag_lim)] + filters = z_filters + Mag_filters + [duty_cycle_key] if self.duty_cycle_on else z_filters + Mag_filters + + # retreive data from mock catalog + catalog_data = catalog_instance.get_quantities([mag_field_key, mag_agn_key], filters=filters) + d = GCRQuery(*((np.isfinite, col) for col in catalog_data)).filter(catalog_data) + + #select point-like sources according to specified agn_flux_fraction + N_all = len(d[mag_field_key]) + fluxid = 'AGN+galaxy' + if self.agn_flux_fraction > 0. and self.agn_flux_fraction < 1.: + point_like_mask = (d[mag_agn_key] - d[mag_field_key] < -2.5*np.log10(self.agn_flux_fraction)) + mags = d[mag_field_key][point_like_mask] + elif self.agn_flux_fraction == 0.: + mags = d[mag_field_key] + else: + mags = d[mag_agn_key] + fluxid = 'AGN' + + ##################################################### + # caclulate the number densities of AGN + ##################################################### + + # get the total number of AGN passing cut and save for txt file + N_tot = len(mags) + if self.agn_flux_fraction > 0.: + total_txt = '{}/{} with AGN magnitude fraction > {}'.format(N_tot, N_all, self.agn_flux_fraction) + fraction_txt = '$\\rm F_{{AGN}}/F_{{Total}} > {}$'.format(self.agn_flux_fraction) + else: + total_txt = '{} AGN'.format(N_tot) + fraction_txt = '' + fraction_txt = '$\\rm F_{{AGN}}/F_{{Total}} > {}$'.format(self.agn_flux_fraction) + + # define the apparent magnitude bins for plotting purposes + dmag = self.validation_data.get('bin_width', 0.25) + + mag_bins = np.arange(self.mag_lo, self.mag_hi + dmag, dmag) + mag_cen = (mag_bins[:-1] + mag_bins[1:])/2 + + # calculate differential binned data (validation data does not divide by bin width) + Ndm, _ = np.histogram(mags, bins=mag_bins) + dn = Ndm/sky_area + + # calculate N( 0) + x = val_mags[mask] + y = np.log10(val_data[mask]) + f = interp1d(x, y, fill_value='extrapolate') + val_mask = (mag_pts > validation_range[0]) & (mag_pts < validation_range[1]) + interp_data = 10**f(mag_pts[val_mask]) + + return mag_pts[val_mask], (cat_data[val_mask] - interp_data)/interp_data + + + def decorate_plot(self, ax, ylabel, scale='log', xlabel=None, legend_title='', legend_loc='upper left'): + """ + """ + ax.set_ylabel(ylabel, size=self.font_size) + leg = ax.legend(loc=legend_loc, fancybox=True, framealpha=0.5, + fontsize=self.legend_size, numpoints=1) + leg.set_title(legend_title, prop={'size':self.legend_size}) + ax.set_yscale(scale) + if xlabel: + ax.set_xlabel(xlabel, size=self.font_size) + else: + for axlabel in ax.get_xticklabels(): + axlabel.set_visible(False) + + + @staticmethod + def post_process_plot(fig): + """ + """ + fig.tight_layout() + fig.subplots_adjust(hspace=0) + fig.subplots_adjust(top=0.94) + + @staticmethod + def save_quantities(keyname, results, filename, comment=''): + """ + """ + if keyname in results: + if keyname+'-' in results and keyname+'+' in results: + fields = ('mag', keyname, keyname+'-', keyname+'+') + header = ', '.join(('Data columns are: ', keyname, keyname+'-', keyname+'+')) + elif keyname+'+-' in results: + fields = ('mag', keyname, keyname+'+-') + header = ', '.join(('Data columns are: ', keyname, keyname+'+-')) + else: + fields = ('mag', keyname) + header = ', '.join(('Data columns are: ', keyname)) + np.savetxt(filename, np.vstack((results[k] for k in fields)).T, fmt='%12.4e', header=': '.join([comment, header])) + + + def conclude_test(self, output_dir): + """ + output_dir: output directory + """ + self.post_process_plot(self.summary_fig) + self.summary_fig.savefig(os.path.join(output_dir, 'summary.png')) + print('saved') + plt.close(self.summary_fig) diff --git a/descqa/configs/AGN_NumberDensity_SDSSg.yaml b/descqa/configs/AGN_NumberDensity_SDSSg.yaml new file mode 100644 index 00000000..ff40a80b --- /dev/null +++ b/descqa/configs/AGN_NumberDensity_SDSSg.yaml @@ -0,0 +1,12 @@ +subclass_name: AGN_NumberDensities.AGN_NumberDensity +band: g +rest_frame_band: i +Mag_lim: -22.5 +z_lim: [0.3, 2.2] +observation: 'SDSS' +validation_range: [15.5, 19.] +agn_flux_fraction: 0. +duty_cycle_on: true +no_agn_extinction: true +included_by_default: true +description: 'Plot N and N(