Skip to content

Commit

Permalink
Merge pull request #185 from LSSTDESC/dndmag
Browse files Browse the repository at this point in the history 
Cosmetic changes for nicer figures
  • Loading branch information
@evevkovacs
evevkovacs committed Sep 30, 2019
2 parents 210c53c + 37f0a5a commit 1772458
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 36 deletions.
97 changes: 61 additions & 36 deletions descqa/apparent_mag_func_test.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,9 @@
from __future__ import unicode_literals, absolute_import, division
import os
import re
import numpy as np
from scipy.interpolate import interp1d
from .utils import get_sky_area
from .base import BaseValidationTest, TestResult
from .plotting import plt

Expand Down Expand Up @@ -41,6 +43,15 @@ def __init__(self, band='r', band_lim=(24.0, 27.5), fractional_tol=0.4, observat
"""
# pylint: disable=super-init-not-called

self.truncate_cat_name = kwargs.get('truncate_cat_name', False)
self.title_in_legend = kwargs.get('title_in_legend', False)
self.skip_label_detail = kwargs.get('skip_label_detail', False)
self.font_size = kwargs.get('font_size', 16)
self.legend_size = kwargs.get('legend_size', 10)
self.x_lower_limit = kwargs.get('x_lower_limit', 15)
self.print_title = kwargs.get('print_title', False)
self.min_mag = kwargs.get('min_mag', 19.)

# catalog quantities needed
possible_mag_fields = ('mag_{}_lsst',
'mag_true_{}_lsst',
Expand All @@ -49,7 +60,8 @@ def __init__(self, band='r', band_lim=(24.0, 27.5), fractional_tol=0.4, observat
'mag_{}_des',
'mag_true_{}_des',
'mag_{}_hsc',
'mag_true_{}_hsc')
'mag_true_{}_hsc',
'mag_{}')
self.possible_mag_fields = [f.format(band) for f in possible_mag_fields]

# attach some attributes to the test
Expand All @@ -58,13 +70,13 @@ def __init__(self, band='r', band_lim=(24.0, 27.5), fractional_tol=0.4, observat
self.fractional_tol = fractional_tol

# set color of lines in plots
colors = plt.cm.jet(np.linspace(0,1,5)) # pylint: disable=no-member
colors = plt.cm.jet(np.linspace(0, 1, 5)) # pylint: disable=no-member
if band == 'g': self.line_color = colors[0]
elif band == 'r': self.line_color = colors[1]
elif band == 'i': self.line_color = colors[2]
elif band == 'z': self.line_color = colors[3]
elif band == 'y': self.line_color = colors[4]
else: self.line_color='black'
else: self.line_color = 'black'

# check for validation observation
if not observation:
Expand All @@ -76,8 +88,8 @@ def __init__(self, band='r', band_lim=(24.0, 27.5), fractional_tol=0.4, observat

# prepare summary plot
self.summary_fig = plt.figure()
upper_rect = 0.2,0.4,0.7,0.55
lower_rect = 0.2,0.125,0.7,0.275
upper_rect = 0.2, 0.4, 0.7, 0.55
lower_rect = 0.2, 0.125, 0.7, 0.275
self.summary_upper_ax, self.summary_lower_ax = self.summary_fig.add_axes(upper_rect), self.summary_fig.add_axes(lower_rect)

def get_validation_data(self, band, observation):
Expand All @@ -96,7 +108,7 @@ def get_validation_data(self, band, observation):
data = np.loadtxt(data_path, unpack=True, usecols=data_args['usecols'], skiprows=data_args['skiprows'])

validation_data = dict(zip(data_args['colnames'], data))
validation_data['label'] = data_args['label']
validation_data['label'] = data_args['label'] if not self.skip_label_detail else data_args['label'].rpartition('(')[0]

return validation_data

Expand All @@ -106,23 +118,27 @@ def post_process_plot(self, upper_ax, lower_ax):
"""

#upper panel
upper_ax.legend(loc='upper left')
upper_ax.set_ylabel(r'$n(< {\rm mag}) ~[{\rm deg^{-2}}]$')
lgnd_title = ''
title = str(self.band_lim[0]) + ' < '+self.band + ' < ' + str(self.band_lim[1])
if self.title_in_legend:
lgnd_title = title
elif self.print_title:
upper_ax.set_title(title)
upper_ax.legend(loc='upper left', title=lgnd_title, fontsize=self.legend_size)
upper_ax.set_ylabel(r'$n(< {\rm mag}) ~[{\rm deg^{-2}}]$', size=self.font_size)
upper_ax.xaxis.set_visible(False)
upper_ax.set_ylim([1000, 10**7])
upper_ax.fill_between([self.band_lim[0], self.band_lim[1]], [0, 0], [10**9, 10**9], alpha=0.1, color='grey')
upper_ax.set_yscale('log')
upper_ax.set_title(str(self.band_lim[0]) + ' < '+self.band + ' < ' + str(self.band_lim[1]))
upper_ax.set_xlim([15,30])
upper_ax.set_xlim([self.x_lower_limit, 30])

#lower panel
lower_ax.fill_between([self.band_lim[0], self.band_lim[1]], [-1, -1], [1, 1], alpha=0.1, color='grey')
lower_ax.set_xlabel(self.band + ' magnitude')
lower_ax.set_ylabel(r'$\Delta n/n$')
lower_ax.set_ylim([-1,1])
lower_ax.set_yticks([-0.6,0.0,0.6])
lower_ax.set_xlim([15,30])

lower_ax.set_xlabel(self.band + ' magnitude', size=self.font_size)
lower_ax.set_ylabel(r'$\Delta n/n$', size=self.font_size)
lower_ax.set_ylim([-1, 1])
lower_ax.set_yticks([-0.6, 0.0, 0.6])
lower_ax.set_xlim([self.x_lower_limit, 30])


def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
Expand All @@ -133,38 +149,45 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
if not mag_field_key:
return TestResult(skipped=True, summary='Catalog is missing requested quantity: {}'.format(self.possible_mag_fields))

# check to see if catalog is a light cone
# this is required since we must be able to calculate the angular area
# if attribute `lightcone` does not exist, allow the catalog to proceed
if not getattr(catalog_instance, 'lightcone', True):
return TestResult(skipped=True, summary="Catalog is not a light cone.")

# obtain or calculate sky area
sky_area = getattr(catalog_instance, 'sky_area', None)
if sky_area is None:
if not catalog_instance.has_quantities(['ra', 'dec']):
return TestResult(skipped=True, summary="'ra' and/or 'dec' not available to compute sky area")
sky_area = get_sky_area(catalog_instance) # compute area from ra and dec

sky_area_label = ' (Sky Area = {:.1f} $\\rm deg^2$)'.format(sky_area)

#####################################################
# caclulate the cumulative number density of galaxies
#####################################################

# filter on extended sources if quantity is available in catalog (eg. in object catalog)
filters = ['extendedness == 1'] if catalog_instance.has_quantity('extendedness') else None

# retreive data from mock catalog
d = catalog_instance.get_quantities([mag_field_key])
d = catalog_instance.get_quantities([mag_field_key], filters=filters)
m = d[mag_field_key]
m = np.sort(m) # put into order--bright to faint

# check to see if catalog is a light cone
# this is required since we must be able to calculate the angular area
if not catalog_instance.lightcone:
return TestResult(skipped=True, summary="Catalog is not a light cone.")

# check to see the angular area if an attribute of the catalog
try:
sky_area = catalog_instance.sky_area
except AttributeError:
return TestResult(skipped=True, summary="Catalog needs an attribute 'sky_area'.")

# get the total number of galaxies in catalog
N_tot = len(m)
N = np.cumsum(np.ones(N_tot))/sky_area

# define the apparent magnitude bins for plotting purposes
dmag = 0.1 # bin widths
max_mag = self.band_lim[1] + 1.0 # go one mag beyond the limit
min_mag = self.band_lim[0] - 1.0 # start at bright galaxies
min_mag = self.min_mag # start at bright galaxies
mag_bins = np.arange(min_mag, max_mag, dmag)

# calculate N(<mag) at the specified points
inds = np.searchsorted(m,mag_bins)
inds = np.searchsorted(m, mag_bins)
mask = (inds >= len(m))
inds[mask] = -1 # take care of edge case
sampled_N = N[inds]
Expand All @@ -174,13 +197,15 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
#################################################

fig = plt.figure()
upper_rect = 0.2,0.4,0.7,0.55
lower_rect = 0.2,0.125,0.7,0.275
upper_rect = 0.2, 0.4, 0.7, 0.55
lower_rect = 0.2, 0.125, 0.7, 0.275
upper_ax, lower_ax = fig.add_axes(upper_rect), fig.add_axes(lower_rect)

# plot on both this plot and any summary plots
upper_ax.plot(mag_bins, sampled_N, '-', label=catalog_name)
self.summary_upper_ax.plot(mag_bins, sampled_N, '-', label=catalog_name)
if self.truncate_cat_name:
catalog_name = re.split('_', catalog_name)[0]
upper_ax.plot(mag_bins, sampled_N, '-', label=catalog_name + sky_area_label)
self.summary_upper_ax.plot(mag_bins, sampled_N, '-', label=catalog_name + sky_area_label)

# plot validation data
n = self.validation_data['n(<mag)']
Expand All @@ -193,7 +218,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
#################################

# interpolate the validation data in order to compare to the mock catalog at same points
non_zero_mask = (self.validation_data['n(<mag)']>0.0)
non_zero_mask = (self.validation_data['n(<mag)'] > 0.0)
x = self.validation_data['mag'][non_zero_mask]
y = np.log10(self.validation_data['n(<mag)'])[non_zero_mask]
f_xy = interp1d(x, y, fill_value='extrapolate')
Expand All @@ -214,7 +239,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
self.summary_lower_ax.plot(mag_bins, delta, '-', label=catalog_name)

# apply 'passing' criterion
if max_frac_diff>self.fractional_tol:
if max_frac_diff > self.fractional_tol:
score = max_frac_diff
passed = False
else:
Expand Down
23 changes: 23 additions & 0 deletions descqa/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,29 @@ def get_sky_volume(sky_area, zlo, zhi, cosmology):
sky_area_rad = np.deg2rad(np.deg2rad(sky_area))
return (dhi**3.0 - dlo**3.0) * sky_area_rad / 3.0

def get_sky_area(catalog_instance, nside=1024):
"""
Parameters
----------
catalog_instance: GCRCatalogs intance
nside: nside parameter for healpy
Returns
-------
sky_area : float
in units of deg**2.0
"""
possible_area_qs = (('ra_true', 'ra'), ('dec_true', 'dec'))
area_qs = [catalog_instance.first_available(*a) for a in possible_area_qs]

pixels = set()
for d in catalog_instance.get_quantities(area_qs, return_iterator=True):
pixels.update(hp.ang2pix(nside, d[area_qs[0]], d[area_qs[1]], lonlat=True))

frac = len(pixels) / hp.nside2npix(nside)
sky_area = frac * np.rad2deg(np.rad2deg(4.0*np.pi))

return sky_area


def get_opt_binpoints(N, sumM, sumM2, bins):
"""
Expand Down

0 comments on commit 1772458

Please sign in to comment.