From d0cf07da031f829132f0a3ec418026ae23caf293 Mon Sep 17 00:00:00 2001 From: Sidney Mau Date: Mon, 10 Jun 2019 14:31:26 -0500 Subject: [PATCH] Migrate from .fits to .npy as per #9 --- README.md | 4 +- config.yaml | 80 +----- simple/diagnostic_plots.py | 419 ------------------------------ simple/farm_plots.py | 3 +- simple/farm_simple.py | 8 +- simple/make_npy.py | 40 +++ simple/make_plot.py | 107 -------- simple/plotting/farm_plots.py | 3 +- simple/plotting/farm_web_plots.py | 3 +- simple/simple_utils.py | 27 +- simple/utils/check_fracdet.py | 3 +- simple/utils/make_webpage.py | 3 +- simple/utils/plot_candidates.py | 3 +- simple/utils/plot_distribution.py | 3 +- simple/utils/sim_check.py | 3 +- 15 files changed, 83 insertions(+), 626 deletions(-) delete mode 100644 simple/diagnostic_plots.py create mode 100644 simple/make_npy.py delete mode 100644 simple/make_plot.py diff --git a/README.md b/README.md index 15dda26..8a7e1f3 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,7 @@ If you have fracdet or maglim files, those should also be copied or linked here; To run the simple binning search, run `farm_simple.py` in `simple_run/`. This will run `search_algorithm.py` over the given data set and write the output to `results_dir/`, logging each job in `results_dir/log_dir`. -The results can then be compiled into a candidate list by running `make_list.py` from `simple_run/` (this is saved as a `.fits` file). +The results can then be compiled into a candidate list by running `make_list.py` from `simple_run/` (this is saved as a `.npy` file). To produce plots from a candidate list, run `farm_plots.py` in `simple_run/`. The output will be written to `save_dir/` with logs in `save_dir/log_dir/`. @@ -43,5 +43,3 @@ The output will be written to `save_dir/` with logs in `save_dir/log_dir/`. By default, `farm_plots.py` will only produce plots for hotspots with statistical significance greater than 5.5 sigma. This threshold has intentionally chosen to be low (such that investigations can be made of very low-significance hotsposts and candidates) but also to minimize junk. - -This code is currently going through an architectural overhaul to capitalize on the OOP nature of python. diff --git a/config.yaml b/config.yaml index b9c5942..40e3af3 100644 --- a/config.yaml +++ b/config.yaml @@ -1,6 +1,5 @@ setup: simple_dir : '/home/s1/smau/software/simple/simple' - #simple_dir : '/Users/keithbechtol/Documents/DES/code/simple' output: results_dir : 'results_dir' @@ -38,52 +37,10 @@ y3_gold_2_0: mag_err: 'SOF_PSF_MAG_ERR_{}' mag_dered: 'MAG_DERED_{}' - mag_g : 'MAG_G' - mag_r : 'MAG_R' - mag_g_err : 'SOF_PSF_MAG_ERR_G' - mag_r_err : 'SOF_PSF_MAG_ERR_R' - #star_galaxy_classification : 'EXTENDED_CLASS_MASH_SOF' - #maglim_g : '/home/s1/kadrlica/projects/y3q2/v6/release/depth/y2q1_maglim_g_n1024_ring.fits.gz' - #maglim_r : '/home/s1/kadrlica/projects/y2q/v6/release/depth/y2q1_maglim_r_n1024_ring.fits.gz' - fracdet_gz : 'y3a2_griz_o.4096_t.32768_coverfoot_EQU.fits.gz' fracdet : 'y3a2_griz_o.4096_t.32768_coverfoot_EQU_decompressed.fits' - -y3a2: - mode : 1 # 0: real, 1: real+sim - sim_population : 'v3/sim_population_v3.fits' - sim_dir: 'v3' - nside : 32 - datadir : '/home/s1/kadrlica/projects/y3a2/dsphs/v2/skim' - candidate_list : 'candidate_list.fits' - mag_max : 24.5 - #candidate_list : '/home/s1/kadrlica/projects/y3a2/dsphs/v2/v2.1/search/ugali_candidates.fits' # ugali - isoname: 'Bressan2012' - isosurvey: 'des' - - basis_1: 'RA' - basis_2: 'DEC' - mag_g : 'MAG_G' - mag_r : 'MAG_R' - mag_g_err : 'PSF_MAG_ERR_G' - mag_r_err : 'PSF_MAG_ERR_R' - - mag_g_dred_flag : 'PSF_MAG_SFD_G' - mag_r_dred_flag : 'PSF_MAG_SFD_R' - mag_err_g_flag : 'PSF_MAG_ERR_G' - mag_err_r_flag : 'PSF_MAG_ERR_R' - flags_g : 'SEXTRACTOR_FLAGS_G' - flags_r : 'SEXTRACTOR_FLAGS_R' - - star_galaxy_classification : 'EXTENDED_CLASS_MASH' - - maglim_g : '/home/s1/kadrlica/projects/y3q2/v6/release/depth/y2q1_maglim_g_n1024_ring.fits.gz' - maglim_r : '/home/s1/kadrlica/projects/y2q/v6/release/depth/y2q1_maglim_r_n1024_ring.fits.gz' - - fracdet : 'y3a2_griz_o.4096_t.32768_coverfoot_EQU.fits.gz' - bliss: # Need to update this mode: 0 sim_population: None @@ -106,13 +63,6 @@ bliss: # Need to update this mag_err: 'WAVG_MAGERR_PSF_{}' mag_dered: 'MAG_DERED_{}' - mag_g : 'MAG_G' - mag_r : 'MAG_R' - mag_g_err : 'WAVG_MAGERR_PSF_G' - mag_r_err : 'WAVG_MAGERR_PSF_R' - #mag_g_err : 'PSF_MAG_ERR_G' - #mag_r_err : 'PSF_MAG_ERR_R' - fracdet : None maglites: # Update @@ -125,22 +75,13 @@ maglites: # Update basis_1: 'RA' basis_2: 'DEC' - mag_g : 'MAG_G' - mag_r : 'MAG_R' - #mag_g_dred_flag : 'PSF_MAG_SFD_G' - #mag_r_dred_flag : 'PSF_MAG_SFD_R' - mag_g_flag : 'WAVG_MAG_PSF_G' - mag_r_flag : 'WAVG_MAG_PSF_R' - mag_err_g_flag : 'WAVG_MAGERR_PSF_G' - mag_err_r_flag : 'WAVG_MAGERR_PSF_R' - extinction_g_flag : 'EXTINCTION_G' - extinction_r_flag : 'EXTINCTION_R' - #star_galaxy_classification : 'EXTENDED_CLASS_MASH' # adapt this for dif data sets - - spread_model_r_flag : 'WAVG_SPREAD_MODEL_R' # use for BISS/MagLiteS - spread_model_r_err_flag : 'SPREADERR_MODEL_R' - flags_g : 'SEXTRACTOR_FLAGS_G' - flags_r : 'SEXTRACTOR_FLAGS_R' + basis_1: 'RA' + basis_2: 'DEC' + band_1: 'G' + band_2: 'R' + mag: 'MAG_PSF_SFD_{}' + mag_err: 'WAVG_MAGERR_PSF_{}' + mag_dered: 'MAG_DERED_{}' fracdet : None @@ -163,12 +104,6 @@ panstarrs: basis_1: 'RA' basis_2: 'DEC' - mag_g : 'MAG_G' - mag_r : 'MAG_R' - mag_g_err : 'GFPSFMAGERR' - mag_r_err : 'RFPSFMAGERR' - flags_g : 'GFLAGS' - flags_r : 'RFLAGS' band_1: 'G' band_2: 'R' @@ -178,4 +113,3 @@ panstarrs: fracdet_gz : 'panstarrs_pseudo_fracdet.fits.gz' fracdet : 'panstarrs_pseudo_fracdet.fits' - #fracdet : '/Users/keithbechtol/Documents/DES/external_data/pan-starrs/simple/sid/simple_v2/panstarrs_pseudo_fracdet.fits.gz' diff --git a/simple/diagnostic_plots.py b/simple/diagnostic_plots.py deleted file mode 100644 index 0d962c2..0000000 --- a/simple/diagnostic_plots.py +++ /dev/null @@ -1,419 +0,0 @@ -#!/usr/bin/env python -""" -Diagnostic plot functions -""" -__author__ = "Sidney Mau" - -import os -import glob -import yaml - -import fitsio as fits -from astropy.coordinates import SkyCoord -from ugali.utils import healpix -from ugali.isochrone import factory as isochrone_factory -import ugali.utils.projector -import ugali.utils.plotting -import healpy - -import numpy as np -from operator import add -from scipy import interpolate -from scipy.signal import argrelextrema -import scipy.ndimage - -import pylab as plt -import matplotlib -from matplotlib import mlab -from mpl_toolkits.axes_grid1 import make_axes_locatable - -import simple.filters -import simple.simple_utils - -################################################################################ - -with open('config.yaml', 'r') as ymlfile: - cfg = yaml.load(ymlfile) - - survey = cfg['survey'] - nside = cfg[survey]['nside'] - datadir = cfg[survey]['datadir'] - isoname = cfg[survey]['isoname'] - isosurvey = cfg[survey]['isosurvey'] - mag_max = cfg[survey]['mag_max'] - basis_1 = cfg[survey]['basis_1'] - basis_2 = cfg[survey]['basis_2'] - - mode = cfg[survey]['mode'] - sim_population = cfg[survey]['sim_population'] - - mag_g = cfg[survey]['mag_g'] - mag_r = cfg[survey]['mag_r'] - mag_g_err = cfg[survey]['mag_g_err'] - mag_r_err = cfg[survey]['mag_r_err'] - - band_1 = cfg[survey]['band_1'] - band_2 = cfg[survey]['band_2'] - mag = cfg[survey]['mag'] - mag_err = cfg[survey]['mag_err'] - mag_dered = cfg[survey]['mag_dered'] - -# construct mags -mag_1 = mag.format(band_1.upper()) -mag_2 = mag.format(band_2.upper()) -mag_err_1 = mag_err.format(band_1.upper()) -mag_err_2 = mag_err.format(band_2.upper()) -mag_dered_1 = mag_dered.format(band_1.upper()) -mag_dered_2 = mag_dered.format(band_2.upper()) - -################################################################################ - -def analysis(targ_ra, targ_dec, mod, mc_source_id): - """Analyze a candidate""" - - pix_nside_select = ugali.utils.healpix.angToPix(nside, targ_ra, targ_dec) - ra_select, dec_select = ugali.utils.healpix.pixToAng(nside, pix_nside_select) - pix_nside_neighbors = np.concatenate([[pix_nside_select], healpy.get_all_neighbours(nside, pix_nside_select)]) - - # Construct data - #data = simple.simple_utils.construct_modal_data(mode, pix_nside_neighbors) - data = simple.simple_utils.construct_real_data(pix_nside_neighbors) - if (mode == 0): - print('mode = 0: running only on real data') - elif (mode == 1): - print('mode = 1: running on real data and simulated data') - - # inject objects for simulated object of mc_source_id - sim_data = simple.simple_utils.construct_sim_data(pix_nside_neighbors, mc_source_id) - data = simple.simple_utils.inject_sim(data, sim_data, mc_source_id) - else: - print('No mode specified; running only on real data') - - print('Loading data...') - data = simple.simple_utils.construct_modal_data(mode, pix_nside_neighbors, mc_source_id) - quality_cut = simple.filters.quality_filter(survey, data) - data = data[quality_cut] - print('Found {} objects...').format(len(data)) - - data = simple.filters.dered_mag(survey, data) - - # This should be generalized to also take the survey - iso = isochrone_factory(name=isoname, survey=isosurvey, age=12, z=0.0001, distance_modulus=mod, band_1=band_1.lower(), band_2=band_2.lower()) - - # g_radius estimate - filter = simple.filters.star_filter(survey, data) - - iso_filter = simple.simple_utils.cut_isochrone_path(data[mag_dered_1], data[mag_dered_2], data[mag_err_1], data[mag_err_2], iso, radius=0.1, return_all=False) - - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data[basis_1], data[basis_2]) - - bins = np.linspace(0, 0.4, 21) # deg - centers = 0.5*(bins[1:] + bins[0:-1]) - area = np.pi*(bins[1:]**2 - bins[0:-1]**2) * 60**2 - hist = np.histogram(angsep[(angsep < 0.4) & filter & iso_filter], bins=bins)[0] # counts - - f_interp = interpolate.interp1d(np.linspace(centers[0], centers[-1], len(hist)), hist/area, 'cubic') - f_range = np.linspace(centers[0], centers[-1], 1000) - f_val = f_interp(f_range) - - pairs = zip(f_range, f_val) - - peak = max(pairs[:len(pairs)/4], key=lambda x: x[1]) # find peak within first quarter - - def peak_index(pairs, peak): - for i in range(len(pairs)): - if pairs[i] == peak: - return i - - osc = int(0.04/0.4*1000) # +/- 0.04 (rounded down) deg oscillation about local extremum - relmin = argrelextrema(f_val, np.less, order=osc)[0] - - try: - if len(relmin) > 0: - #half_point = f_range[relmin[0]] - i = 0 - while ((f_range[relmin[i]] <= f_range[peak_index(pairs,peak)]) & (i <= len(relmin)-1)): - i+=1 - half_point = f_range[relmin[i]] - elif len(relmin) == 0: - half_peak = (peak[1] + np.mean(f_val[len(f_val)/4:]))/2. # normalized to background (after first quarter) - #half_peak = np.mean(f_val[len(f_val)/4:]) - half_pairs = [] - for i in pairs[peak_index(pairs, peak):len(pairs)/2]: # start after peak, stay within first quarter - if i != peak: - half_pairs.append((i[0], abs(i[1]-half_peak))) - half_point = min(half_pairs, key=lambda x: x[1])[0] # deg - except: - half_point = 0.1 # fixed value to catch errors - - g_min = 0.5/60. # deg - g_max = 12./60. # deg - - if half_point < g_min: - g_radius = g_min - elif half_point > g_max: - g_radius = g_max - else: - g_radius = half_point # deg - - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data[basis_1], data[basis_2]) - nbhd = (angsep < g_radius) - - return(data, iso, g_radius, nbhd) - -def density_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, type): - """Stellar density plot""" - - if type == 'stars': - filter = simple.filters.star_filter(survey, data) - plt.title('Stellar Density') - elif type == 'galaxies': - filter = simple.filters.galaxy_filter(survey, data) - plt.title('Galactic Density') - elif type == 'blue_stars': - filter = simple.filters.color_filter(survey, data) \ - & simple.filters.star_filter(survey, data) - plt.title('Blue Stellar Density') - - iso_filter = simple.simple_utils.cut_isochrone_path(data[mag_dered_1], data[mag_dered_2], data[mag_err_1], data[mag_err_2], iso, radius=0.1, return_all=False) - - # projection of image - proj = ugali.utils.projector.Projector(targ_ra, targ_dec) - x, y = proj.sphereToImage(data[filter & iso_filter][basis_1], data[filter & iso_filter][basis_2]) - - bound = 0.5 #1. - steps = 100. - bins = np.linspace(-bound, bound, steps) - - signal = np.histogram2d(x, y, bins=[bins, bins])[0] - - sigma = 0.01 * (0.25 * np.arctan(0.25*g_radius*60. - 1.5) + 1.3) - - convolution = scipy.ndimage.filters.gaussian_filter(signal, sigma/(bound/steps)) - plt.pcolormesh(bins, bins, convolution.T, cmap='Greys') - - plt.xlim(bound, -bound) - plt.ylim(-bound, bound) - plt.gca().set_aspect('equal') - plt.xlabel(r'$\Delta \alpha$ (deg)') - plt.ylabel(r'$\Delta \delta$ (deg)') - - ax = plt.gca() - divider = make_axes_locatable(ax) - cax = divider.append_axes('right', size = '5%', pad=0) - plt.colorbar(cax=cax) - -def star_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd): - """Star bin plot""" - - filter = simple.filters.star_filter(survey, data) - - iso_filter = simple.simple_utils.cut_isochrone_path(data[mag_dered_1], data[mag_dered_2], data[mag_err_1], data[mag_err_2], iso, radius=0.1, return_all=False) - - # projection of image - proj = ugali.utils.projector.Projector(targ_ra, targ_dec) - x, y = proj.sphereToImage(data[filter & iso_filter][basis_1], data[filter & iso_filter][basis_2]) - - plt.scatter(x, y, edgecolor='none', s=3, c='black') - plt.xlim(0.2, -0.2) - plt.ylim(-0.2, 0.2) - plt.gca().set_aspect('equal') - plt.xlabel(r'$\Delta \alpha$ (deg)') - plt.ylabel(r'$\Delta \delta$ (deg)') - - plt.title('Stars') - -def cm_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, type): - """Color-magnitude plot""" - - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data[basis_1], data[basis_2]) - annulus = (angsep > g_radius) & (angsep < 1.) - - if type == 'stars': - filter = simple.filters.star_filter(survey, data) - plt.title('Stellar Color-Magnitude') - elif type == 'galaxies': - filter = simple.filters.galaxy_filter(survey, data) - plt.title('Galactic Color-Magnitude') - - iso_filter = simple.simple_utils.cut_isochrone_path(data[mag_dered_1], data[mag_dered_2], data[mag_err_1], data[mag_err_2], iso, radius=0.1, return_all=False) - - # Plot background objects - plt.scatter(data[mag_dered_1][filter & annulus] - data[mag_dered_2][filter & annulus], data[mag_dered_1][filter & annulus], c='k', alpha=0.1, edgecolor='none', s=1) - - # Plot isochrone - ugali.utils.plotting.drawIsochrone(iso, lw=2, label='{} Gyr, z = {}'.format(iso.age, iso.metallicity)) - - # Plot objects in nbhd - plt.scatter(data[mag_dered_1][filter & nbhd] - data[mag_dered_2][filter & nbhd], data[mag_dered_2][filter & nbhd], c='g', s=5, label='r < {:.3f}$^\circ$'.format(g_radius)) - - # Plot objects in nbhd and near isochrone - plt.scatter(data[mag_dered_1][filter & nbhd & iso_filter] - data[mag_dered_2][filter & nbhd & iso_filter], data[mag_dered_1][filter & nbhd & iso_filter], c='r', s=5, label='$\Delta$CM < 0.1') - - plt.axis([-0.5, 1, 16, mag_max]) - plt.gca().invert_yaxis() - plt.gca().set_aspect(1./4.) - plt.legend(loc='upper left') - plt.xlabel('{} - {} (mag)'.format(band_1.lower(), band_2.lower())) - plt.ylabel('{} (mag)'.format(band_1.lower())) - -def hess_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd): - """Hess plot""" - - filter = simple.filters.star_filter(survey, data) - - plt.title('Hess') - - c1 = SkyCoord(targ_ra, targ_dec, unit='deg') - - r_near = 2.*g_radius # annulus begins at 2*g_radius away from centroid - r_far = np.sqrt(5.)*g_radius # annulus has same area as inner area - - #inner = (c1.separation(SkyCoord(data[basis_1], data[basis_2], unit='deg')).deg < g_radius) - #outer = (c1.separation(SkyCoord(data[basis_1], data[basis_2], unit='deg')).deg > r_near) & (c1.separation(SkyCoord(data[basis_1], data[basis_2], unit='deg')).deg < r_far) - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data[basis_1], data[basis_2]) - inner = (angsep < g_radius) - outer = ((angsep > r_near) & (angsep < r_far)) - - xbins = np.arange(-0.5, 1.1, 0.1) - ybins = np.arange(16., mag_max + 0.5, 0.5) - - foreground = np.histogram2d(data[mag_dered_1][inner & filter] - data[mag_dered_2][inner & filter], data[mag_dered_1][inner & filter], bins=[xbins, ybins]) - background = np.histogram2d(data[mag_dered_1][outer & filter] - data[mag_dered_2][outer & filter], data[mag_dered_1][outer & filter], bins=[xbins, ybins]) - - fg = foreground[0].T - bg = background[0].T - - fg_abs = np.absolute(fg) - bg_abs = np.absolute(bg) - - mask_abs = fg_abs + bg_abs - mask_abs[mask_abs == 0.] = np.nan # mask signficiant zeroes - - signal = fg - bg - signal = np.ma.array(signal, mask=np.isnan(mask_abs)) # mask nan - - cmap = matplotlib.cm.viridis - cmap.set_bad('w', 1.) - plt.pcolormesh(xbins, ybins, signal, cmap=cmap) - - plt.colorbar() - - ugali.utils.plotting.drawIsochrone(iso, lw=2, c='k', zorder=10, label='Isocrhone') - - plt.axis([-0.5, 1.0, 16, mag_max]) - plt.gca().invert_yaxis() - plt.gca().set_aspect(1./4.) - plt.xlabel('{} - {} (mag)'.format(band_1.lower(), band_2.lower())) - plt.ylabel('{} (mag)'.format(band_1.lower())) - -def radial_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, field_density=None): - """Radial distribution plot""" - - ## Deprecated? - #filter_s = simple.filters.star_filter(survey, data) - #filter_g = simple.filters.galaxy_filter(survey, data) - - plt.title('Radial Distribution') - - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data[basis_1], data[basis_2]) - - # Isochrone filtered/unfiltered - iso_seln_f = simple.simple_utils.cut_isochrone_path(data[mag_dered_1], data[mag_dered_2], data[mag_err_1], data[mag_err_2], iso, radius=0.1, return_all=False) - iso_seln_u = ~iso_seln_f - - bins = np.linspace(0, 0.4, 21) # deg - centers = 0.5*(bins[1:] + bins[0:-1]) - area = np.pi*(bins[1:]**2 - bins[0:-1]**2) * 60**2 - - def interp_values(type, seln): - if type == 'stars': - filter = simple.filters.star_filter(survey, data) - elif type == 'galaxies': - filter = simple.filters.galaxy_filter(survey, data) - - if seln == 'f': - iso_filter = iso_seln_f - elif seln == 'u': - iso_filter = iso_seln_u - - hist = np.histogram(angsep[(angsep < 0.4) & filter & iso_filter], bins=bins)[0] # counts - - f_interp = interpolate.interp1d(np.linspace(centers[0], centers[-1], len(hist)), hist/area, 'cubic') - f_range = np.linspace(centers[0], centers[-1], 1000) - f_val = f_interp(f_range) - - return(f_range, f_val) - - def value_errors(type, seln): - if type == 'stars': - filter = simple.filters.star_filter(survey, data) - elif type == 'galaxies': - filter = simple.filters.galaxy_filter(survey, data) - if seln == 'f': - iso_filter = iso_seln_f - elif seln == 'u': - iso_filter = iso_seln_u - - hist = np.histogram(angsep[(angsep < 0.4) & filter & iso_filter], bins=bins)[0] # counts - - val = hist/area - yerr = np.sqrt(hist)/area - - return(val, yerr) - - f_range, f_val = interp_values('stars', 'f') - pairs = zip(f_range, f_val) - peak = max(pairs[:len(pairs)/4], key=lambda x: x[1]) # find peak within first quarter - def peak_index(pairs, peak): - for i in range(len(pairs)): - if pairs[i] == peak: - return i - - plt.axvline(x=f_range[peak_index(pairs,peak)], color='m', label='peak') - - plt.axvline(x=g_radius, color='r', label='g_radius') - - f_range, f_val = interp_values('galaxies', 'f') - plt.plot(f_range, f_val, '-g', label='Filtered Galaxies') - - f_range, f_val = interp_values('stars', 'u') - plt.plot(f_range, f_val, '-k', alpha=0.25, label='Unfiltered Stars') - - val, y_err = value_errors('stars', 'f') - plt.plot(centers, val, '.b') - plt.errorbar(centers, val, yerr=y_err, fmt='none', ecolor='b', elinewidth=1, capsize=5) - - f_range, f_val = interp_values('stars', 'f') - plt.plot(f_range, f_val, '-b', label='Filtered Stars') - - if field_density: - plt.axhline(y=field_density, color='blue', ls='--', label='Model Field') - - ymax = plt.ylim()[1] - plt.annotate(r'$\approx %0.1f$' + str(round(g_radius, 3)) + '$^\circ$', (g_radius*1.1, ymax/50.), color='red', bbox=dict(boxstyle='round,pad=0.0', fc='white', alpha=0.75, ec='white', lw=0)) - plt.xlim(bins[0], bins[-1]) - plt.ylim(0., ymax) - plt.legend(loc='upper right') - plt.xlabel('Angular Separation (deg)') - plt.ylabel('Denisty (arcmin$^{-2})$') - -#def maglim_plot(targ_ra, targ_dec, data, iso, band): -# """Maglim plots""" -# -# reso = 0.5 -# xsize = 2.*60./reso -# -# if band == 'g': -# reader = pyfits.open(maglim_g) -# m_maglim_g = reader[1].data.field('I').flatten() -# reader.close() -# m_maglim_g[np.isnan(m_maglim_g)] = healpy.UNSEEN -# #healpy.gnomview(m_maglim_g, fig='summary', rot=(targ_ra, targ_dec, 0.), reso=reso, xsize=xsize, title='maglim g (S/N =10)', sub=(3, 4, 11)) -# healpy.gnomview(m_maglim_g, rot=(targ_ra, targ_dec, 0.), reso=reso, xsize=xsize, title='maglim g (S/N =10)', sub=(3, 4, 8)) -# elif band == 'r': -# reader = pyfits.open(maglim_r) -# m_maglim_r = reader[1].data.field('I').flatten() -# reader.close() -# m_maglim_r[np.isnan(m_maglim_r)] = healpy.UNSEEN -# healpy.gnomview(m_maglim_r, rot=(targ_ra, targ_dec, 0.), reso=reso, xsize=xsize, title='maglim r (S/N =10)', sub=(3, 4, 12)) diff --git a/simple/farm_plots.py b/simple/farm_plots.py index d785c37..0e55ff9 100644 --- a/simple/farm_plots.py +++ b/simple/farm_plots.py @@ -45,7 +45,8 @@ print('Plotting hotspots with sig > {}'.format(sig_cut)) -candidate_list = fits.read(candidate_list) +#candidate_list = fits.read(candidate_list) +candidate_list = np.load(candidate_list) try: # simple candidate_list = candidate_list[candidate_list['SIG'] > sig_cut] except: # ugali diff --git a/simple/farm_simple.py b/simple/farm_simple.py index 0ad5ad8..a6462f1 100644 --- a/simple/farm_simple.py +++ b/simple/farm_simple.py @@ -49,17 +49,17 @@ def submit_job(ra, dec, pix, mc_source_id, mode, **population_file): if (mode == 0): - outfile = '{}/results_nside_{}_{}.csv'.format(results_dir, nside, pix) + outfile = '{}/results_nside_{}_{}.npy'.format(results_dir, nside, pix) logfile = '{}/results_nside_{}_{}.log'.format(log_dir, nside, pix) #command = 'python {}/search_algorithm.py {:0.2f} {:0.2f} {:0.2f} {} {}'.format(simple_dir, ra, dec, mc_source_id, outfile, logfile) elif (mode == 1): - #outfile = '{}/results_mc_source_id_{}.csv'.format(results_dir, mc_source_id) # all values in mc_source_id_array should be the same + #outfile = '{}/results_mc_source_id_{}.npy'.format(results_dir, mc_source_id) # all values in mc_source_id_array should be the same #logfile = '{}/results_mc_source_id_{}.log'.format(log_dir, mc_source_id) # all values in mc_source_id_array should be the same - outfile = '{}/results_nside_{}_{}.csv'.format(results_dir, nside, pix) + outfile = '{}/results_nside_{}_{}.npy'.format(results_dir, nside, pix) logfile = '{}/results_nside_{}_{}.log'.format(log_dir, nside, pix) #command = 'python {}/search_algorithm.py {:0.2f} {:0.2f} {:0.2f} {} {} {}'.format(simple_dir, ra, dec, mc_source_id, outfile, logfile, population_file) elif (mode == 2): - outfile = '{}/results_mc_source_id_{}.csv'.format(results_dir, mc_source_id) # all values in mc_source_id_array should be the same + outfile = '{}/results_mc_source_id_{}.npy'.format(results_dir, mc_source_id) # all values in mc_source_id_array should be the same logfile = '{}/results_mc_source_id_{}.log'.format(log_dir, mc_source_id) # all values in mc_source_id_array should be the same #batch = 'csub -n {} -o {} '.format(jobs, logfile) batch = 'csub -n {} -o {} --host all '.format(jobs, logfile) # testing condor updates diff --git a/simple/make_npy.py b/simple/make_npy.py new file mode 100644 index 0000000..978469e --- /dev/null +++ b/simple/make_npy.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python +""" +Compile candidate list from results_dir +""" +__author__ = "Sidney Mau" + +import glob +import yaml + +import numpy as np + +with open('config.yaml', 'r') as ymlfile: + cfg = yaml.load(ymlfile) + + survey = cfg['survey'] + + basis_1 = cfg[survey]['basis_1'] + basis_2 = cfg[survey]['basis_2'] + candidate_list = cfg[survey]['candidate_list'] + + +fs = glob.glob('{}/*.npy'.format(cfg['output']['results_dir'])) +data = np.concatenate([np.load(f) for f in fs]) +np.save(candidate_list, data) + +# Diagnostic output +data = np.load(candidate_list) +print("{} hotspots found.").format(len(data)) +cut_0 = (data['SIG'] > 5.5) +print("{} hotspots found with SIG > 5.5.").format(len(data[cut_0])) +cut_1 = (data['SIG'] > 10) +print("{} hotspots found with SIG > 10.").format(len(data[cut_1])) +cut_2 = (data['SIG'] > 15) +print("{} hotspots found with SIG > 15.").format(len(data[cut_2])) +cut_3 = (data['SIG'] > 20) +print("{} hotspots found with SIG > 20.").format(len(data[cut_3])) +cut_4 = (data['SIG'] > 25) +print("{} hotspots found with SIG > 25.").format(len(data[cut_4])) +cut_5 = (data['SIG'] >= 37.5) +print("{} hotspots found with SIG >= 37.55").format(len(data[cut_5])) diff --git a/simple/make_plot.py b/simple/make_plot.py deleted file mode 100644 index 5aaa7eb..0000000 --- a/simple/make_plot.py +++ /dev/null @@ -1,107 +0,0 @@ -#!/usr/bin/env python -""" -Arrange and produce plots -""" -__author__ = "Sidney Mau" - -# Set the backend first! -import matplotlib -matplotlib.use('Agg') - -import sys -import os -import yaml -import numpy as np -import numpy -import pylab as plt -from matplotlib import gridspec - -import ugali.utils.projector -import ugali.candidate.associate - -import simple.diagnostic_plots - -print(matplotlib.get_backend()) - -with open('config.yaml', 'r') as ymlfile: - cfg = yaml.load(ymlfile) - - survey = cfg['survey'] - nside = cfg[survey]['nside'] - datadir = cfg[survey]['datadir'] - basis_1 = cfg[survey]['basis_1'] - basis_2 = cfg[survey]['basis_2'] - -save_dir = os.path.join(os.getcwd(), cfg['output']['save_dir']) -if not os.path.exists(save_dir): - os.mkdir(save_dir) - -################################################################# - -try: - targ_ra, targ_dec, mod, sig, mc_source_id, field_density, = float(sys.argv[1]), float(sys.argv[2]), float(sys.argv[3]), float(sys.argv[4]), float(sys.argv[5]), float(sys.argv[6]) -except: - try: - targ_ra, targ_dec, mod, sig, mc_source_id, = float(sys.argv[1]), float(sys.argv[2]), float(sys.argv[3]), float(sys.argv[4]), float(sys.argv[5]) - field_density = None - except: - sys.exit('ERROR! Coordinates not given in correct format.') - -fig = plt.figure(figsize=(20, 17)) -fig.subplots_adjust(wspace=0.3, hspace=0.3) -gs = gridspec.GridSpec(3, 3) - -data, iso, g_radius, nbhd = simple.diagnostic_plots.analysis(targ_ra, targ_dec, mod, mc_source_id) - -print('Making diagnostic plots for ({}, {}) = ({}, {})...'.format(basis_1, basis_2, targ_ra, targ_dec)) - -fig.add_subplot(gs[0,0]) -simple.diagnostic_plots.density_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'stars') - -fig.add_subplot(gs[1,0]) -simple.diagnostic_plots.density_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'galaxies') - -fig.add_subplot(gs[2,0]) -simple.diagnostic_plots.density_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'blue_stars') - -fig.add_subplot(gs[0,1]) -simple.diagnostic_plots.cm_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'stars') - -fig.add_subplot(gs[1,1]) -simple.diagnostic_plots.cm_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'galaxies') - -fig.add_subplot(gs[0,2]) -simple.diagnostic_plots.hess_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd) - -fig.add_subplot(gs[1,2]) -simple.diagnostic_plots.star_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd) - -fig.add_subplot(gs[2,1:3]) -simple.diagnostic_plots.radial_plot(targ_ra, targ_dec, data, iso, g_radius, nbhd, field_density) - -# Name -try: # ugali - association_string = candidate_list[candidate]['NAME'] -except: # simple - # Check for possible associations - glon_peak, glat_peak = ugali.utils.projector.celToGal(targ_ra, targ_dec) - catalog_array = ['McConnachie15', 'Harris96', 'Corwen04', 'Nilson73', 'Webbink85', 'Kharchenko13', 'WEBDA14','ExtraDwarfs','ExtraClusters'] - catalog = ugali.candidate.associate.SourceCatalog() - for catalog_name in catalog_array: - catalog += ugali.candidate.associate.catalogFactory(catalog_name) - - idx1, idx2, sep = catalog.match(glon_peak, glat_peak, tol=0.5, nnearest=1) - match = catalog[idx2] - if len(match) > 0: - association_string = '{} at {:0.3f} deg'.format(match[0]['name'], float(sep)) - else: - association_string = 'No association within 0.5 deg' - -association_string = str(np.char.strip(association_string)) -association_string = repr(association_string) - -plt.suptitle('{}\n'.format(association_string) + r'($\alpha$, $\delta$, $\mu$, $\sigma$, MC_SOURCE_ID) = ({:0.2f}, {:0.2f}, {:0.2f}, {:0.2f}, {:0.0f})'.format(targ_ra, targ_dec, mod, sig, mc_source_id), fontsize=24) - -file_name = 'candidate_{:0.2f}_{:0.2f}'.format(targ_ra, targ_dec) -plt.savefig(save_dir+'/'+file_name+'.png', bbox_inches='tight') -plt.close() diff --git a/simple/plotting/farm_plots.py b/simple/plotting/farm_plots.py index 3f78b27..d66a8e6 100644 --- a/simple/plotting/farm_plots.py +++ b/simple/plotting/farm_plots.py @@ -45,7 +45,8 @@ print('Plotting hotspots with sig > {}'.format(sig_cut)) -candidate_list = fits.read(candidate_list) +#candidate_list = fits.read(candidate_list) +candidate_list = np.load(candidate_list) try: # simple candidate_list = candidate_list[candidate_list['SIG'] > sig_cut] except: # ugali diff --git a/simple/plotting/farm_web_plots.py b/simple/plotting/farm_web_plots.py index f22eb0d..bf2c810 100644 --- a/simple/plotting/farm_web_plots.py +++ b/simple/plotting/farm_web_plots.py @@ -45,7 +45,8 @@ print('Plotting hotspots with sig > {}'.format(sig_cut)) -candidate_list = fits.read(candidate_list) +#candidate_list = fits.read(candidate_list) +candidate_list = np.load(candidate_list) try: # simple candidate_list = candidate_list[candidate_list['SIG'] > sig_cut] except: # ugali diff --git a/simple/simple_utils.py b/simple/simple_utils.py index 41e7eba..9512dbe 100644 --- a/simple/simple_utils.py +++ b/simple/simple_utils.py @@ -696,18 +696,21 @@ def search_by_object(nside, data, distance_modulus, pix_nside_select, ra_select, def write_output(results_dir, nside, pix_nside_select, ra_peak_array, dec_peak_array, r_peak_array, distance_modulus_array, n_obs_peak_array, n_obs_half_peak_array, n_model_peak_array, sig_peak_array, mc_source_id_array, mode, outfile): - writer = open(outfile, 'a') # append if exists - for ii in range(0, len(sig_peak_array)): - # SIG, RA, DEC, MODULUS, r, n_obs, n_model, mc_source_id - writer.write('{:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}\n'.format(sig_peak_array[ii], - ra_peak_array[ii], - dec_peak_array[ii], - distance_modulus_array[ii], - r_peak_array[ii], - n_obs_peak_array[ii], - n_obs_half_peak_array[ii], - n_model_peak_array[ii], - mc_source_id_array[ii])) + #writer = open(outfile, 'a') # append if exists + #for ii in range(0, len(sig_peak_array)): + # # SIG, RA, DEC, MODULUS, r, n_obs, n_model, mc_source_id + # writer.write('{:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}, {:10.3f}\n'.format(sig_peak_array[ii], + # ra_peak_array[ii], + # dec_peak_array[ii], + # distance_modulus_array[ii], + # r_peak_array[ii], + # n_obs_peak_array[ii], + # n_obs_half_peak_array[ii], + # n_model_peak_array[ii], + # mc_source_id_array[ii])) + data = [tuple(row) for row in np.stack([sig_peak_array, ra_peak_array, dec_peak_array, distance_modulus_array, r_peak_array, n_obs_peak_array, n_obs_half_peak_array, n_model_peak_array, mc_source_id_array], axis=-1)] + arr = np.array(data, dtype=[('SIG', float), ('RA', float), ('DEC', float), ('DISTANCE_MODULUS', float), ('R_PEAK', float), ('N_OBS_PEAK', float), ('N_OBS_HALF_PEAK', float), ('N_MODEl_PEAK', float), ('MC_SOURCE_ID', float)]) + np.save(outfile, arr) ######################################################################## diff --git a/simple/utils/check_fracdet.py b/simple/utils/check_fracdet.py index d73b6af..89a9ed9 100644 --- a/simple/utils/check_fracdet.py +++ b/simple/utils/check_fracdet.py @@ -24,7 +24,8 @@ candidate_list = cfg[survey]['candidate_list'] -data = fits.read(candidate_list) +#data = fits.read(candidate_list) +data = np.load(candidate_list) data = data[data['SIG'] > 15] #fracdet = fits.read('{}_pseudo_fracdet.fits.gz'.format(survey)) diff --git a/simple/utils/make_webpage.py b/simple/utils/make_webpage.py index b1d38b5..9788f7b 100644 --- a/simple/utils/make_webpage.py +++ b/simple/utils/make_webpage.py @@ -54,7 +54,8 @@ #candidate_list = np.genfromtxt('remains_des_simple.txt', names=True) #candidate_list[::-1].sort(order='SIG') -candidate_list = fits.read(candidate_list) +#candidate_list = fits.read(candidate_list) +candidate_list = np.load(candidate_list) candidate_list[::-1].sort(order='SIG') ################################################################################ diff --git a/simple/utils/plot_candidates.py b/simple/utils/plot_candidates.py index d492e29..a372183 100644 --- a/simple/utils/plot_candidates.py +++ b/simple/utils/plot_candidates.py @@ -26,7 +26,8 @@ if not os.path.exists(log_dir): os.mkdir(log_dir) -data = fits.read(candidate_list) +#data = fits.read(candidate_list) +data = np.load(candidate_list) try: sig_sel, plot = float(sys.argv[1]), float(sys.argv[2]) diff --git a/simple/utils/plot_distribution.py b/simple/utils/plot_distribution.py index 0047876..15b6188 100644 --- a/simple/utils/plot_distribution.py +++ b/simple/utils/plot_distribution.py @@ -26,7 +26,8 @@ basis_1 = cfg[survey]['basis_1'] basis_2 = cfg[survey]['basis_2'] -data = fits.read(candidate_list) +#data = fits.read(candidate_list) +data = np.load(candidate_list) #data = data[data['SIG'] > 15] plt.figure() diff --git a/simple/utils/sim_check.py b/simple/utils/sim_check.py index 77a22ef..b23f097 100644 --- a/simple/utils/sim_check.py +++ b/simple/utils/sim_check.py @@ -60,7 +60,8 @@ ######################################################################## -data = fits.read(candidate_list) +#data = fits.read(candidate_list) +data = np.load(candidate_list) sim_pop = fits.read(sim_population) #data = data[(sim_pop['DIFFICULTY'] == 0) & (sim_pop['N_CATALOG'] > 0)]