diff --git a/ugali/simple/README.md b/ugali/simple/README.md deleted file mode 100644 index 697e800..0000000 --- a/ugali/simple/README.md +++ /dev/null @@ -1,9 +0,0 @@ -# Simple Binning Search Program - -This code has been adapted from Keith Bechtol's original simple binning code program to be more modular and ready-to-use for future searches. - -[Sidney Mau](https://github.com/SidneyMau) - -## Configuration and use - -The `config.yaml` file will point the different parts of the code to the data set over which the search will be performed and the directories where the results will be saved. Once this has been set, then `farm_simple.py` will use `search_algorithm` to perform the simple binnning search over the given data set, storing the results in `results_dir/`. Then, `make_list.py` will compile the data stored in `results_dir/` into a candidate list, `results.csv`. Finally, running `render_plots.py` will use `diagnostic_plots.py` to create diagnostic plots for each candidate with sigma > 5.5 and will save these plots in `save_dir/`. diff --git a/ugali/simple/config.yaml b/ugali/simple/config.yaml deleted file mode 100644 index d6867d3..0000000 --- a/ugali/simple/config.yaml +++ /dev/null @@ -1,36 +0,0 @@ -#data: 'y3a2' - -#y3a2: -nside : 32 -#datadir : '/home/s1/kadrlica/projects/y3a2/dsphs/v1/skim' -datadir : '/home/s1/kadrlica/projects/y3a2/dsphs/v2/skim' -results_dir : 'results_dir' -log_dir : 'log_dir' -save_dir : 'save_dir_y3a2_low_sig' -candidate_list : 'candidate_list.csv' - -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 : None - -mag_g_dred_flag : 'PSF_MAG_SFD_G' -mag_r_dred_flag : 'PSF_MAG_SFD_R' -#mag_g_flag : 'CM_MAG_G' -#mag_r_flag : 'CM_MAG_R' -mag_err_g_flag : 'PSF_MAG_ERR_G' -mag_err_r_flag : 'PSF_MAG_ERR_R' -#extinction_g_flag : 'EXTINCTION_G' -#extinction_r_flag : 'EXTINCTION_R' -star_galaxy_classification : 'EXTENDED_CLASS_MASH' -#star_filter : '' -#spread_model_r_flag : 'CM_T' # use for BISS/MagLiteS -#spread_model_r_err_flag : 'CM_T_ERR' -flags_g : 'SEXTRACTOR_FLAGS_G' -flags_r : 'SEXTRACTOR_FLAGS_R' - -# star, galaxy, blue_star filters defined with a string of selections -# make another config string that takes in kind of data (i.e. DES, BLISS, MagLiteS) and then determines which selections to put together in which order to make the appropriate filter -# also include sub flags for reddened vs de-redenned to handle how those strings should work -# this will put a fair bit of emphasis on the structure of the config file, but it may help generalize the code to these various data sets. Then, any new format of data (i.e. LSST) can just used a new "class" within the config file. -# Ultimately, this should separate the data handling from the analysis of the data diff --git a/ugali/simple/diagnostic_plots.py b/ugali/simple/diagnostic_plots.py deleted file mode 100644 index f4c355c..0000000 --- a/ugali/simple/diagnostic_plots.py +++ /dev/null @@ -1,466 +0,0 @@ -#!/usr/bin/env python -""" -Diagnostic plot functions -""" -__author__ = "Sidney Mau" - -import os - -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 utils - -import pylab as plt -import astropy.io.fits as pyfits -import matplotlib -from matplotlib import mlab -from mpl_toolkits.axes_grid1 import make_axes_locatable - -import yaml - -################################################################################ - -with open('config.yaml', 'r') as ymlfile: - cfg = yaml.load(ymlfile) - -nside = cfg['nside'] -datadir = cfg['datadir'] - -maglim_g = cfg['maglim_g'] -maglim_r = cfg['maglim_r'] - -mag_g_dred_flag = cfg['mag_g_dred_flag'] -mag_r_dred_flag = cfg['mag_r_dred_flag'] -#mag_g_flag = cfg['mag_g_flag'] -#mag_r_flag = cfg['mag_r_flag'] -mag_err_g_flag = cfg['mag_err_g_flag'] -mag_err_r_flag = cfg['mag_err_r_flag'] -#extinction_g_flag = cfg['extinction_g_flag'] -#extinction_r_flag = cfg['extinction_r_flag'] -star_galaxy_classification = cfg['star_galaxy_classification'] -#spread_model_r_flag = cfg['spread_model_r_flag'] -#spread_model_r_err_flag = cfg['spread_model_r_err_flag'] -flags_g = cfg['flags_g'] -flags_r = cfg['flags_r'] - -#COLUMNS = ['RA', 'DEC'] -#COLUMNS+= ['WAVG_MAG_PSF_G', 'WAVG_MAG_PSF_R'] -#COLUMNS+= ['EXTINCTION_G', 'EXTINCTION_R'] -#COLUMNS+= ['WAVG_SPREAD_MODEL_R', 'SPREADERR_MODEL_R'] -#COLUMNS+= ['WAVG_MAGERR_PSF_G', 'WAVG_MAGERR_PSF_R'] -#COLUMNS+= ['MAGERR_PSF_G', 'MAGERR_PSF_R'] - -################################################################################ - -def star_filter(data): - """Selects stars from the data set""" - filter = (data[star_galaxy_classification] >= 0) & (data[star_galaxy_classification] <= 1) - #filter = (data[spread_model_r_flag] < 0.003) - return filter - -def galaxy_filter(data): - """Selects galaxies from the data set""" - #filter = (data[spread_model_r_flag] > 0.005) - filter = (data[star_galaxy_classification] > 1) - return filter - -def blue_star_filter(data): - """Selects blue stars from the data set""" - filter = (star_filter(data)) & ((data[mag_g_dred_flag] - data[mag_r_dred_flag]) < 0.4) # 0.2 - return filter - -#filter_r = mag_r <= 23 # Filters out background galaxy contamination -#filter not working propery - -################################################################################ - -def analysis(targ_ra, targ_dec, mod): - """Analyze a candidate""" - - #hdisc = healpix.ang2disc(nside, targ_ra, targ_dec, 1, inclusive=True) - #files = [] - #for i in hdisc: - # if os.path.exists(datadir+'/cat_hpx_'+str(i).zfill(5)+'.fits'): - # files.append(datadir+'/cat_hpx_'+str(i).zfill(5)+'.fits') - # - ##data = utils.load_infiles(files, columns=['RA', 'DEC', 'WAVG_MAG_PSF_G', 'WAVG_MAG_PSF_R', 'EXTINCTION_G', 'EXTINCTION_R', 'WAVG_SPREAD_MODEL_R', 'SPREADERR_MODEL_R', 'WAVG_MAGERR_PSF_G', 'WAVG_MAGERR_PSF_R', 'MAGERR_PSF_G', 'MAGERR_PSF_R']) - #data = utils.load_infiles(files, columns=COLUMNS) - - 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)]) - data_array = [] - for pix_nside in pix_nside_neighbors: - #infile = '%s/cat_hpx_%05i.fits'%(datadir, pix_nside) - infile = '%s/y3a2_ngmix_cm_%05i.fits'%(datadir, pix_nside) - #infile = '%s/*_%05i.fits'%(datadir, pix_nside) # TODO: get to work - if not os.path.exists(infile): - continue - reader = pyfits.open(infile) - data_array.append(reader[1].data) - reader.close() - print('Assembling data...') - data = np.concatenate(data_array) - print('Found {} objects...').format(len(data)) - print('Loading data...') - - ## De-redden magnitudes - #try: - # data = mlab.rec_append_fields(data, ['WAVG_MAG_PSF_DRED_G', 'WAVG_MAG_PSF_DRED_R'], [data[mag_g_dred_flag], data[mag_r_dred_flag]]) - ##except: - ## data = mlab.rec_append_fields(data, ['WAVG_MAG_PSF_DRED_G', 'WAVG_MAG_PSF_DRED_R'], [data[mag_g_flag] - data[extinction_g_flag], data[mag_r_flag] - data[extinction_r_flag]]) - #except: - # data = mlab.rec_append_fields(data, ['WAVG_MAG_PSF_DRED_G', 'WAVG_MAG_PSF_DRED_R'], [data[mag_g_flag], data[mag_r_flag]]) - # - #mag_g = data['WAVG_MAG_PSF_DRED_G'] - #mag_r = data['WAVG_MAG_PSF_DRED_R'] - mag_g = data[mag_g_dred_flag] - mag_r = data[mag_r_dred_flag] - - iso = isochrone_factory('Bressan2012', age=12, z=0.0001, distance_modulus=mod) - - # g_radius estimate - filter_s = star_filter(data) - - mag_g = data[mag_g_dred_flag] - mag_r = data[mag_r_dred_flag] - - iso_filter = (iso.separation(mag_g, mag_r) < 0.1) - - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data['RA'], data['DEC']) - - 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_s & 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 = list(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]] # TODO rename - 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 - - #c1 = SkyCoord(targ_ra, targ_dec, unit='deg') # frame is ICRS - #nbhd = c1.separation(SkyCoord(data['RA'], data['DEC'], unit='deg')).deg < g_radius # selects objects inside the galactic radius - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data['RA'], data['DEC']) - nbhd = (angsep < g_radius) - - return(data, iso, g_radius, nbhd) - -def densityPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, type): - """Stellar density plot""" - - mag_g = data[mag_g_dred_flag] - mag_r = data[mag_r_dred_flag] - - if type == 'stars': - filter = star_filter(data) - plt.title('Stellar Density') - elif type == 'galaxies': - filter = galaxy_filter(data) - plt.title('Galactic Density') - elif type == 'blue_stars': - filter = blue_star_filter(data) - plt.title('Blue Stellar Density') - - iso_filter = (iso.separation(mag_g, mag_r) < 0.1) - - # projection of image - proj = ugali.utils.projector.Projector(targ_ra, targ_dec) - x, y = proj.sphereToImage(data[filter & iso_filter]['RA'], data[filter & iso_filter]['DEC']) # filter & iso_filter - - 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) # full range, arctan - - 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 starPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd): - """Star bin plot""" - - mag_g = data[mag_g_dred_flag] - mag_r = data[mag_r_dred_flag] - - filter = star_filter(data) - - iso_filter = (iso.separation(mag_g, mag_r) < 0.1) - - # projection of image - proj = ugali.utils.projector.Projector(targ_ra, targ_dec) - x, y = proj.sphereToImage(data[filter & iso_filter]['RA'], data[filter & iso_filter]['DEC']) - - 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 cmPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, type): - """Color-magnitude plot""" - - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data['RA'], data['DEC']) - annulus = (angsep > g_radius) & (angsep < 1.) - - mag_g = data[mag_g_dred_flag] - mag_r = data[mag_r_dred_flag] - - if type == 'stars': - filter = star_filter(data) - plt.title('Stellar Color-Magnitude') - elif type == 'galaxies': - filter = galaxy_filter(data) - plt.title('Galactic Color-Magnitude') - - iso_filter = (iso.separation(mag_g, mag_r) < 0.1) - - # Plot background objects - plt.scatter(mag_g[filter & annulus] - mag_r[filter & annulus], mag_g[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(mag_g[filter & nbhd] - mag_r[filter & nbhd], mag_g[filter & nbhd], c='g', s=5, label='r < {:.3f}$^\circ$'.format(g_radius)) - - # Plot objects in nbhd and near isochrone - plt.scatter(mag_g[filter & nbhd & iso_filter] - mag_r[filter & nbhd & iso_filter], mag_g[filter & nbhd & iso_filter], c='r', s=5, label='$\Delta$CM < 0.1') - - plt.axis([-0.5, 1, 16, 24]) - plt.gca().invert_yaxis() - plt.gca().set_aspect(1./4.) - plt.legend(loc='upper left') - plt.xlabel('g-r (mag)') - plt.ylabel('g (mag)') - -def hessPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd): - """Hess plot""" - - mag_g = data[mag_g_dred_flag] - mag_r = data[mag_r_dred_flag] - - filter_s = star_filter(data) - - plt.title('Hess') - - c1 = SkyCoord(targ_ra, targ_dec, unit='deg') - - r_near = 2.*g_radius # annulus begins at 3*g_radius away from centroid - r_far = np.sqrt(5.)*g_radius # annulus has same area as inner area - - inner = (c1.separation(SkyCoord(data['RA'], data['DEC'], unit='deg')).deg < g_radius) - outer = (c1.separation(SkyCoord(data['RA'], data['DEC'], unit='deg')).deg > r_near) & (c1.separation(SkyCoord(data['RA'], data['DEC'], unit='deg')).deg < r_far) - - xbins = np.arange(-0.5, 1.1, 0.1) - ybins = np.arange(16., 24.5, 0.5) - - foreground = np.histogram2d(mag_g[inner & filter_s] - mag_r[inner & filter_s], mag_g[inner & filter_s], bins=[xbins, ybins]) - background = np.histogram2d(mag_g[outer & filter_s] - mag_r[outer & filter_s], mag_g[outer & filter_s], 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, 24]) - plt.gca().invert_yaxis() - plt.gca().set_aspect(1./4.) - plt.xlabel('g-r (mag)') - plt.ylabel('g (mag)') - - #ax = plt.gca() - #divider = make_axes_locatable(ax) - #cax = divider.append_axes('right', size = '5%', pad=0) - #plt.colorbar(cax=cax) - -def radialPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd): - """Radial distribution plot""" - - mag_g = data[mag_g_dred_flag] - mag_r = data[mag_r_dred_flag] - - filter_s = star_filter(data) - filter_g = galaxy_filter(data) - - plt.title('Radial Distribution') - - angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data['RA'], data['DEC']) - - # Isochrone filtered/unfiltered - iso_seln_f = (iso.separation(mag_g, mag_r) < 0.1) - iso_seln_u = (iso.separation(mag_g, mag_r) >= 0.1) - - 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 = star_filter(data) - elif type == 'galaxies': - filter = galaxy_filter(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 = star_filter(data) - elif type == 'galaxies': - filter = galaxy_filter(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 = list(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') - - 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/ugali/simple/extra/compile_results_v3.py b/ugali/simple/extra/compile_results_v3.py deleted file mode 100644 index 437a466..0000000 --- a/ugali/simple/extra/compile_results_v3.py +++ /dev/null @@ -1,318 +0,0 @@ -import glob -import sys -import os -import numpy -import pylab -import healpy -import pyfits - -import ugali.utils.healpix -#import uagli.utils.projector -import ugali.candidate.associate - -pylab.ion() - -############################################################ -""" -datadir = '/Users/keithbechtol/Documents/DES/projects/calibration/Y2N/data/catalog' -infiles = glob.glob('%s/starcat_y2n_gr_v3_00000*.fits'%(datadir)) -data_array = [] -for infile in infiles: - print infile - reader = pyfits.open(infile) - data_array.append(reader[1].data) - reader.close() -print 'Assembling data...' -data = numpy.concatenate(data_array) - -############################################################ - -print 'Applying cuts...' -cut = (data['WAVG_MAG_PSF_G'] < 23.5) \ - & ((data['WAVG_MAG_PSF_G'] - data['WAVG_MAG_PSF_R']) < 1.) \ - & (numpy.fabs(data['SPREAD_MODEL_R']) < 0.003 + data['SPREADERR_MODEL_R']) -#& (numpy.fabs(data['WAVG_SPREAD_MODEL_R']) < 0.003 + data['SPREADERR_MODEL_R']) -data = data[cut] - -############################################################ - -print 'Pixelizing...' -npix_256 = healpy.nside2npix(256) -pix_256 = ugali.utils.healpix.angToPix(256, data['RA'], data['DEC']) -m_256 = numpy.histogram(pix_256, numpy.arange(npix_256 + 1))[0].astype(float) -m_256[m_256 == 0.] = healpy.UNSEEN -""" -############################################################ - -save = False - -infiles = glob.glob('results_midway_v8/*.csv') - -ra_array = [] -dec_array = [] -sig_array = [] -r_array = [] -distance_modulus_array = [] -for infile in infiles: - data = numpy.recfromcsv(infile) - if data.shape == (0,): - continue - - if data['ra'].shape == (): - ra_array.append([data['ra']]) - dec_array.append([data['dec']]) - sig_array.append([data['sig']]) - r_array.append([data['r']]) - distance_modulus_array.append([data['distance_modulus']]) - else: - ra_array.append(data['ra']) - dec_array.append(data['dec']) - sig_array.append(data['sig']) - r_array.append(data['r']) - distance_modulus_array.append(data['distance_modulus']) - -sig_array = numpy.concatenate(sig_array) -index_sort = numpy.argsort(sig_array) -sig_array = sig_array[index_sort] -ra_array = numpy.concatenate(ra_array)[index_sort] -dec_array = numpy.concatenate(dec_array)[index_sort] -r_array = numpy.concatenate(r_array)[index_sort] -distance_modulus_array = numpy.concatenate(distance_modulus_array)[index_sort] - -### Run Associations ### - -catalog_array = ['McConnachie12', '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) -association_array = [] -association_angsep_array = numpy.tile(180., len(sig_array)) -for ii in range(0, len(sig_array)): - glon, glat = ugali.utils.projector.celToGal(ra_array[ii], dec_array[ii]) - idx1, idx2, angsep = catalog.match(glon, glat, tol=0.5, nnearest=1) - match = catalog[idx2] - if len(match) > 0: - association_array.append(match[0]['name'].replace(' ', '_')) - association_angsep_array[ii] = angsep - else: - association_array.append('NONE') -association_array = numpy.array(association_array) - - -#cut = ((r_array < 0.21) | (r_array > 0.31)) & (sig_array > 5.9) & (distance_modulus_array < 23.)# & (distance_modulus_array > 20.) -#cut = (r_array < 0.2) & (sig_array > 7.) & (dec_array < -40.) -#cut = (r_array < 0.2) & (sig_array > 5.6) & (sig_array < 5.7) -#cut = (r_array > 0.2) & (sig_array > 10) -#cut = (r_array > 0.2) & (r_array < 0.28) & (sig_array > 9.)# & (sig_array < 7.) -#cut = (r_array < 0.2) & (sig_array > 7.) -#cut = (r_array < 0.2) & (sig_array > 6.) & (sig_array < 7.) -#cut = (r_array < 0.2) & (sig_array > 5.5) & (sig_array < 6.) -cut = numpy.logical_not((r_array > 0.2) | \ - (((numpy.char.count(association_array, 'NGC_') > 0) | \ - (numpy.char.count(association_array, 'UGC_') > 0) | \ - (numpy.char.count(association_array, 'IC_') > 0)) & (association_angsep_array < 0.2))) & (sig_array > 5.5) - -#cut = ugali.utils.projector.angsep(ra_array, dec_array, 55.20, -37.50) < 0.1 -#cut = ugali.utils.projector.angsep(ra_array, dec_array, 56.36, -60.44) < 0.1 -#cut = ugali.utils.projector.angsep(ra_array, dec_array, 67.01, -44.34) < 0.1 - -sig_cut_array = [5.5, 6., 6.5, 7., 10.] -for sig_cut in sig_cut_array: - print(sig_cut, numpy.sum((sig_array > sig_cut) & (r_array < 0.2))) - -for ii in range(0, numpy.sum(cut)): - print(ra_array[cut][ii], dec_array[cut][ii], sig_array[cut][ii], association_array[cut][ii], association_angsep_array[cut][ii]) - -for ra, dec, distance_modulus in zip(ra_array[cut], dec_array[cut], distance_modulus_array[cut]): - outfile = 'candidate_%.2f_%.2f.png'%(ra, dec) - if distance_modulus < 30.: - pass - os.system('cp figs_midway_v8/%s figs_temp/.'%(outfile)) - #os.system('open figs_midway_v1/%s'%(outfile)) - -#sys.exit('DONE') - -""" -pylab.figure() -#pylab.scatter(ra_array[cut], dec_array[cut], c=sig_array[cut]) -pylab.scatter(ra_array, dec_array, c=sig_array, edgecolor='none') -colorbar = pylab.colorbar() -colorbar.set_label(r'Significance ($\sigma$)') -""" - -gc_catalog = ugali.candidate.associate.catalogFactory('Harris96') -sat_catalog = ugali.candidate.associate.catalogFactory('McConnachie12') - -ra_ngc, dec_ngc = list(zip(*[[54.63, -35.45], - [11.78, -20.76], - [51.59, -21.34], - [28.25, -13.74], - [50.59, -37.26], - [51.62, -35.72], - [319.12, -48.36], - [45.40, -14.85], - [66.93, -55.02]])) - -ra_des, dec_des, name_des = list(zip(*[[53.92, -54.05, 'Ret II'], - [56.09, -43.53, 'Eri II'], - [343.06, -58.57, 'Tuc II'], - [43.87, -54.11, 'Hor I'], - [317.20, -51.16, 'Kim 2'], - [70.95, -50.28, 'Pic I'], - [354.99, -54.41, 'Phe II'], - [35.69, -52.28, 'Eri III'], - [344.18, -50.16, 'Gru I'], - [49.13, -50.02, 'Hor II'], - [8.51, -49.04, 'Luque 1']])) - -ra_des_new, dec_des_new, name_des_new = list(zip(*[[331.06, -46.47, 'Gru II'], - [359.04, -59.63, 'Tuc III'], - [82.85, -28.03, 'Col I'], # Rock solid - #[32.29, -12.17, 'Cet II'], - [0.668, -60.888, 'Tuc IV'], - [354.37, -63.26, 'Tuc V'], # Also solid - #[70.95, -50.29, 'New'], # Pic I - #[94.85, -50.33, 'New'], - #[67.01, -44.34, 'New'], - [56.36, -60.44, 'Many faint stars'], - [30.02, 3.35, 'Wide near S82']])) - #[8.51, -49.04, 'Luque 1'], - #[344.19, -50.16, 'New'], # Grus I - #[317.21, -51.16, 'New'], # Ind I - #[55.20, -37.50, 'New'], - #[14.08, -18.36, 'New'], - #[348.73, -42.56, 'New']]) - -# Cannot immediately rule out -#55.20, -37.50 # Looks OK, but there is a somewhat bright neaby galaxy in Aladin viewer, needs further investigation -#317.21, -51.16 # Very compact, very few stars, no association, clean in Aladin viewer, Ind I -#81.18, -54.69 # Not super compelling, but possibly a low surface brightness overdensity -#14.08, -18.36 # Very sparse if there at all, but stars line up nicely along isochrone, overdensity visible in map, no association, clean in Aladin viewer -#21.76, -6.04 # Right on edge of footprint, so a bit hard to tell -#24.37, -34.44 # Nothing obviously wrong, but not many stars at MSTO -#26.83, -51.93 # Spotty -#30.37, -3.26 # Whiting 1 -#344.19, -50.16 # Looks pretty good, no association, clean in Aladin viewer, Grus I -#348.73, -42.56 # A good example of something at a distance that makes very difficult to distinguish from foreground, but an apparent overdensity -#38.80, 3.84 # Sketchy -#43.90, -60.05 # Overdensity is somewhat visible in the map, MSTO just below g ~ 23 -#56.36, -60.44 # Clear overdensity, no association, clean in Aladin viewer -#58.76, -49.61 # AM 1 -#64.70, -24.09 # Sparse, but slight overdensity visible -#66.19, -21.19 # Eri -#67.01, -44.34 # Clearly visible overdensity, very sparse, but stars line up on isochrone, no association, clean in Aladin viewer -#70.95, -50.29 # Looks very good, no association, clean in Aladin viewer, Pic I -#94.85, -50.33 # Possibly very near, no association, clean in Aladin viewer - -#[343.04, -58.59, 'New 2'], Tuc II -#[43.89, -54.12, 'New 3'], Hor I - -#ra_des, dec_des, name_des = zip(*[[359.04, -59.63, 'Tuc III']]) - -#m_256 = numpy.load('density_map_nside_512.npy') -m_256 = numpy.load('/Users/keithbechtol/Documents/DES/projects/mw_substructure/des/y2n/skymap/density_map_nside_1024_v6.npy') -m_256[m_256 == 0.] = healpy.UNSEEN - -pylab.figure(num='map1', figsize=(16, 10)) -healpy.mollview(m_256, - fig='map1', cmap='binary', xsize=6400, min=0.5 * numpy.median(m_256[m_256 > 0.]), max=4. * numpy.median(m_256[m_256 > 0.]), - unit='Stellar Density', - title='Inclusive') -pylab.xlim(-1, 0.55) -pylab.ylim(-0.85, 0.2) -healpy.projscatter(ra_array, dec_array, c=sig_array, edgecolor='none', lonlat=True, vmin=5., vmax=25.) -if save: - pylab.savefig('compile_map_inclusive.png', dpi=150, bbox_inches='tight') - -pylab.figure(num='map2', figsize=(16, 10)) -healpy.mollview(m_256, - fig='map2', cmap='binary', xsize=6400, min=0.5 * numpy.median(m_256[m_256 > 0.]), max=4. * numpy.median(m_256[m_256 > 0.]), - unit='Stellar Density', - title=r'Cleaned: size < 0.2 deg & significance > 6 $\sigma$ & m - M < 23.') -pylab.xlim(-1, 0.55) -pylab.ylim(-0.85, 0.2) -healpy.projscatter(ra_array[cut], dec_array[cut], c=sig_array[cut], edgecolor='none', lonlat=True, vmin=5., vmax=25.) -#healpy.projscatter(ra_array, dec_array, c=sig_array, edgecolor='none', lonlat=True, vmin=5., vmax=25.) -healpy.projscatter(ra_des, dec_des, edgecolor='red', facecolor='none', marker='o', s=40, lonlat=True) -healpy.projscatter(ra_des_new, dec_des_new, edgecolor='blue', facecolor='none', marker='o', s=50, lonlat=True) -healpy.projscatter(gc_catalog.data['ra'], gc_catalog.data['dec'], edgecolor='red', facecolor='none', marker='x', s=40, lonlat=True) -healpy.projscatter(sat_catalog.data['ra'], sat_catalog.data['dec'], edgecolor='blue', facecolor='none', marker='x', s=40, lonlat=True) -healpy.projscatter(ra_ngc, dec_ngc, edgecolor='green', facecolor='none', marker='x', s=40, lonlat=True) -if save: - pylab.savefig('compile_map_clean.png', dpi=150, bbox_inches='tight') - - -#pylab.figure() -#pylab.scatter(ra_array, dec_array, c=sig_array, edgecolor='none', vmin=5., vmax=25.) -#pylab.scatter(ra_array[cut], dec_array[cut], c=sig_array[cut], edgecolor='none', vmin=5., vmax=25.) -#pylab.colorbar() -#pylab.scatter(ra_des, dec_des, c='black', edgecolor='black', marker='x') - - -pylab.figure() -#pylab.hist(distance_modulus_array[(r_array < 0.2) & (sig_array > 8.)], bins=numpy.arange(16.25, 24.25, 0.5)) -pylab.hist(distance_modulus_array, bins=numpy.arange(16.25, 24.25, 0.5)) -pylab.xlabel('m - M (mag)') -pylab.ylabel('N') -pylab.savefig('compile_hist_distance_modulus.png', dpi=150, bbox_inches='tight') - -pylab.figure() -pylab.hist(r_array, bins=numpy.arange(0.005, 0.315, 0.01)) -pylab.xlabel('Size (deg)') -pylab.ylabel('N') -pylab.xlim(0., 0.30) -if save: - pylab.savefig('compile_hist_size.png', dpi=150, bbox_inches='tight') - -cut_dirty = (r_array > 0.2) | \ - (((numpy.char.count(association_array, 'NGC_') > 0) | \ - (numpy.char.count(association_array, 'UGC_') > 0) | \ - (numpy.char.count(association_array, 'IC_') > 0)) & (association_angsep_array < 0.2)) - -pylab.figure() -pylab.yscale('log', nonposy='clip') -pylab.ylim(0.1, 1.e3) -pylab.hist(sig_array, bins=numpy.linspace(5, 25, 41), color='green', alpha=0.5, label='All') -pylab.hist(sig_array[~cut_dirty], bins=numpy.linspace(5, 25, 41), color='green', label='r < 0.2 deg & no large galaxy association') -pylab.xlabel(r'Significance ($\sigma$)') -pylab.ylabel('N') -pylab.legend(loc='upper right', frameon=False) -if save: - pylab.savefig('compile_hist_significance.png', dpi=150, bbox_inches='tight') - #pylab.savefig('figs_summary_v8/compile_hist_significance.png', dpi=150, bbox_inches='tight') - -pylab.figure() -pylab.scatter(sig_array, r_array, c=distance_modulus_array, edgecolor='none') -colorbar = pylab.colorbar() -colorbar.set_label('m - M (mag)') -pylab.xlabel(r'Significance ($\sigma$)') -pylab.ylabel('Size (deg)') -pylab.xlim(4., 26.) -pylab.ylim(0., 0.32) -if save: - pylab.savefig('compile_scatter.png', dpi=150, bbox_inches='tight') - -for ii in range(0, len(ra_des)): - angsep = ugali.utils.projector.angsep(ra_des[ii], dec_des[ii], ra_array, dec_array) - outfile = 'candidate_%.2f_%.2f.png'%(ra_array[numpy.argmin(angsep)], dec_array[numpy.argmin(angsep)]) - #print name_des[ii], sig_array[numpy.argmin(angsep)], outfile - print('| %s | %.2f, %.2f | %.2f | {{thumbnail(%s, size=400)}} |'%(name_des[ii], ra_des[ii], dec_des[ii], sig_array[numpy.argmin(angsep)], outfile)) - os.system('cp figs_midway_v8/%s figs_reported/.'%(outfile)) - -for ii in range(0, len(ra_des_new)): - angsep = ugali.utils.projector.angsep(ra_des_new[ii], dec_des_new[ii], ra_array, dec_array) - outfile = 'candidate_%.2f_%.2f.png'%(ra_array[numpy.argmin(angsep)], dec_array[numpy.argmin(angsep)]) - #print name_des_new[ii], ra_des_new[ii], dec_des_new[ii], sig_array[numpy.argmin(angsep)], outfile - print('| %s | %.2f, %.2f | %.2f | {{thumbnail(%s, size=400)}} |'%(name_des_new[ii], ra_des_new[ii], dec_des_new[ii], sig_array[numpy.argmin(angsep)], outfile)) - os.system('cp figs_midway_v8/%s figs_seed/.'%(outfile)) - - -outfile = 'simple_binner_compiled_v8.csv' -writer = open(outfile, 'w') -writer.write('sig, ra, dec, distance_modulus, r, association, association_angsep\n') -for ii in range(0, len(sig_array)): - writer.write('%10.2f, %10.2f, %10.2f, %10.2f, %10.2f, %30s, %10.2f\n'%(sig_array[ii], - ra_array[ii], dec_array[ii], distance_modulus_array[ii], - r_array[ii], - association_array[ii], association_angsep_array[ii])) -writer.close() diff --git a/ugali/simple/extra/density_map.py b/ugali/simple/extra/density_map.py deleted file mode 100644 index 11d9341..0000000 --- a/ugali/simple/extra/density_map.py +++ /dev/null @@ -1,40 +0,0 @@ -import glob -import numpy -import healpy -import pyfits -import pylab - -import ugali.utils.healpix - -pylab.ion() - -#datadir = '/project/kicp/bechtol/des/mw_substructure/y2n/data/catalog/hpx/cat' -#datadir = '/project/kicp/bechtol/des/mw_substructure/y2n/data/catalog/v4/hpx' -datadir = '/project/kicp/bechtol/des/mw_substructure/y2n/data/catalog/v6/hpx' -infiles = glob.glob('%s/cat_hpx_*.fits'%(datadir)) - -nside = 2**10 -npix = healpy.nside2npix(nside) - -pix_array = [] -for infile in infiles: - print(infile) - reader = pyfits.open(infile) - cut = (reader[1].data['FLAGS_G'] < 4) & (reader[1].data['FLAGS_R'] < 4) \ - & (reader[1].data['QSLR_FLAG_G'] == 0) & (reader[1].data['QSLR_FLAG_R'] == 0) \ - & (reader[1].data['WAVG_MAG_PSF_G'] < 23.) \ - & ((reader[1].data['WAVG_MAG_PSF_G'] - reader[1].data['WAVG_MAG_PSF_R']) < 1.) \ - & (numpy.fabs(reader[1].data['WAVG_SPREAD_MODEL_R']) < 0.003 + reader[1].data['SPREADERR_MODEL_R']) - pix_array.append(ugali.utils.healpix.angToPix(nside, reader[1].data.field('RA')[cut], reader[1].data.field('DEC')[cut])) - reader.close() - -pix_array = numpy.concatenate(pix_array) - -m = numpy.histogram(pix_array, numpy.arange(npix + 1))[0].astype(float) -#m[m == 0] = healpy.UNSEEN - -#numpy.save('density_map_nside_%i_v4.npy'%(nside), m) -numpy.save('density_map_nside_%i_v6.npy'%(nside), m) - -#m = numpy.tile(healpy.UNSEEN) - diff --git a/ugali/simple/farm_simple.py b/ugali/simple/farm_simple.py deleted file mode 100644 index 1e60dc9..0000000 --- a/ugali/simple/farm_simple.py +++ /dev/null @@ -1,58 +0,0 @@ -#!/usr/bin/env python -""" -Perform simple binning search -""" -__author__ = "Sidney Mau" - -import os -import time -import subprocess -import glob -import astropy.io.fits as pyfits -import healpy -import numpy - -import ugali.utils.healpix - -import yaml - -############################################################ - -with open('config.yaml', 'r') as ymlfile: - cfg = yaml.load(ymlfile) - -nside = cfg['nside'] -datadir = cfg['datadir'] - -results_dir = os.path.join(os.getcwd(), cfg['results_dir']) -if not os.path.exists(results_dir): - os.mkdir(results_dir) - -log_dir = os.path.join(os.getcwd(), cfg['results_dir'], cfg['log_dir']) -if not os.path.exists(log_dir): - os.mkdir(log_dir) - -#infiles = glob.glob('%s/cat_hpx_*.fits'%(datadir)) -infiles = glob.glob ('%s/y3a2_ngmix_cm_*.fits'%(datadir)) -############################################################ - -print('Pixelizing...') -pix_nside = [] # Equatorial coordinates, RING ordering scheme -for infile in infiles: - pix_nside.append(int(infile.split('.fits')[0].split('_')[-1])) - -############################################################ - -for ii in range(0, len(pix_nside)): - ra, dec = ugali.utils.healpix.pixToAng(nside, pix_nside[ii]) - - print('({}/{})').format(ii, len(pix_nside)) - - #pix_nside[ii] = pix_nside_select - logfile = '%s/results_nside_%s_%i.log'%(log_dir, nside, pix_nside[ii]) - batch = 'csub -n 20 -o %s '%logfile # q local for debugging - command = 'python search_algorithm.py %.2f %.2f'%(ra, dec) - command_queue = batch + command - print(command_queue) - #os.system('./' + command) # Run locally - os.system(command_queue) # Submit to queue diff --git a/ugali/simple/make_list.py b/ugali/simple/make_list.py deleted file mode 100644 index 2a46cb5..0000000 --- a/ugali/simple/make_list.py +++ /dev/null @@ -1,23 +0,0 @@ -#!/usr/bin/env python -""" -Compile candidate list from results_dir -""" -__author__ = "Sidney Mau" - -import glob -import yaml - -with open('config.yaml', 'r') as ymlfile: - cfg = yaml.load(ymlfile) - -candidate_list = cfg['candidate_list'] - - -results_file = open(candidate_list, 'w') -results_file.write('sig, ra, dec, distance_modulus, r\n') -for file in glob.glob('{}/*.csv'.format(cfg['results_dir'])): - writer = open(file, 'r') - for line in writer: - results_file.write(line) - writer.close() -results_file.close() diff --git a/ugali/simple/render_plots.py b/ugali/simple/render_plots.py deleted file mode 100644 index 6685c92..0000000 --- a/ugali/simple/render_plots.py +++ /dev/null @@ -1,115 +0,0 @@ -#!/usr/bin/env python -""" -Arrange and produce plots -""" -__author__ = "Sidney Mau" - -# Set the backend first! -import matplotlib -matplotlib.use('Agg') - -import os -import yaml -import numpy as np -import numpy -import pylab as plt -from matplotlib import gridspec -from multiprocessing import Pool - -import ugali.utils.projector -import ugali.candidate.associate -import diagnostic_plots - -print(matplotlib.get_backend()) - -with open('config.yaml', 'r') as ymlfile: - cfg = yaml.load(ymlfile) - -nside = cfg['nside'] -datadir = cfg['datadir'] -candidate_list = cfg['candidate_list'] - -save_dir = os.path.join(os.getcwd(), cfg['save_dir']) -if not os.path.exists(save_dir): - os.mkdir(save_dir) - -candidate_list = np.genfromtxt(candidate_list, delimiter=',', names=['sig', 'ra', 'dec', 'distance_modulus', 'r'])[1:] #, 'association', 'association_angsep'])[1:] -candidate_list = candidate_list[candidate_list['sig'] > 5.5] # only plot hotspots of sufficent significance -# try with sig > 10., 5.5, etc... - -################################################################# - -print('{} candidates found...').format(len(candidate_list)) - -def renderPlot(candidate): - """ - Make plots - """ - - fig = plt.figure(figsize=(20, 17)) - fig.subplots_adjust(wspace=0.3, hspace=0.3) - gs = gridspec.GridSpec(3, 3) - - print('Analyzing candidate {}/{}...').format(candidate+1, len(candidate_list)) - - #sig, targ_ra, targ_dec, mod, r = candidate_list[candidate] - sig = round(candidate_list[candidate]['sig'], 2) - targ_ra = round(candidate_list[candidate]['ra'], 2) - targ_dec = round(candidate_list[candidate]['dec'], 2) - mod = candidate_list[candidate]['distance_modulus'] - data, iso, g_radius, nbhd = diagnostic_plots.analysis(targ_ra, targ_dec, mod) - - print('Making diagnostic plots for (RA, Dec) = ({}, {})...').format(targ_ra, targ_dec) - - fig.add_subplot(gs[0,0]) - diagnostic_plots.densityPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'stars') - - fig.add_subplot(gs[1,0]) - diagnostic_plots.densityPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'galaxies') - - fig.add_subplot(gs[2,0]) - diagnostic_plots.densityPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'blue_stars') - - fig.add_subplot(gs[0,1]) - diagnostic_plots.cmPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'stars') - - fig.add_subplot(gs[1,1]) - diagnostic_plots.cmPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, 'galaxies') - - fig.add_subplot(gs[0,2]) - diagnostic_plots.hessPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd) - - fig.add_subplot(gs[1,2]) - diagnostic_plots.starPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd) - - fig.add_subplot(gs[2,1:3]) - diagnostic_plots.radialPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd) - - #try: - # Check for possible associations - glon_peak, glat_peak = ugali.utils.projector.celToGal(targ_ra, targ_dec) - catalog_array = ['McConnachie12', '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 {:.3f} deg'.format(match[0]['name'], float(sep)) - else: - association_string = 'No association within 0.5 deg' - #except: - # association_string = 'Association search error' - - plt.suptitle('{}\n'.format(association_string) + r'($\alpha$, $\delta$, $\mu$, $\sigma$) = ({}, {}, {}, {})'.format(targ_ra, targ_dec, mod, sig), fontsize=24) - - file_name = 'candidate_{}_{}'.format(targ_ra, targ_dec) - plt.savefig(save_dir+'/'+file_name+'.png', bbox_inches='tight') - plt.close() - -#renderPlot(0) -if __name__ == '__main__': - pool = Pool(20) - index = list(range(len(candidate_list))) - pool.map(renderPlot, index) diff --git a/ugali/simple/search_algorithm.py b/ugali/simple/search_algorithm.py deleted file mode 100755 index 48c7aeb..0000000 --- a/ugali/simple/search_algorithm.py +++ /dev/null @@ -1,603 +0,0 @@ -#!/usr/bin/env python -""" -Simple binning search algorithm -""" -__author__ = "Keith Bechtol, Sidney Mau" - -# Set the backend first! -import matplotlib -matplotlib.use('Agg') # May need to disable this option for interactive plotting -import pylab - -import sys -import os -import glob -import numpy as np -from matplotlib import mlab -#import astropy.io.fits as pyfits -import astropy.io.fits as pyfits -import healpy -import scipy.interpolate -from scipy import interpolate -import scipy.ndimage -import scipy.signal -import scipy.stats -import scipy.spatial - -import ugali.utils.healpix -import ugali.utils.projector -import ugali.isochrone -import ugali.utils.plotting -import ugali.candidate.associate - -from ugali.isochrone import factory as isochrone_factory -from astropy.coordinates import SkyCoord - -import yaml - -#pylab.ion() # For interactive plotting - -########################################################### - -#with open('config_dev_v1.yaml', 'r') as ymlfile: -with open('config.yaml', 'r') as ymlfile: - cfg = yaml.load(ymlfile) - -nside = cfg['nside'] -datadir = cfg['datadir'] - -maglim_g = cfg['maglim_g'] -maglim_r = cfg['maglim_r'] - -mag_g_dred_flag = cfg['mag_g_dred_flag'] -mag_r_dred_flag = cfg['mag_r_dred_flag'] -#mag_g_flag = cfg['mag_g_flag'] -#mag_r_flag = cfg['mag_r_flag'] -mag_err_g_flag = cfg['mag_err_g_flag'] -mag_err_r_flag = cfg['mag_err_r_flag'] -#extinction_g_flag = cfg['extinction_g_flag'] -#extinction_r_flag = cfg['extinction_r_flag'] -star_galaxy_classification = cfg['star_galaxy_classification'] -#spread_model_r_flag = cfg['spread_model_r_flag'] -#spread_model_r_err_flag = cfg['spread_model_r_err_flag'] -flags_g = cfg['flags_g'] -flags_r = cfg['flags_r'] - -results_dir = os.path.join(os.getcwd(), cfg['results_dir']) -if not os.path.exists(results_dir): - os.mkdir(results_dir) - -########################################################### - -print(matplotlib.get_backend()) - -############################################################ - -def cutIsochronePath(g, r, g_err, r_err, isochrone, radius=0.1, return_all=False): - """ - Cut to identify objects within isochrone cookie-cutter. - """ - if np.all(isochrone.stage == 'Main'): - # Dotter case - index_transition = len(isochrone.stage) - else: - # Other cases - index_transition = np.nonzero(isochrone.stage > 3)[0][0] + 1 - mag_1_rgb = isochrone.mag_1[0: index_transition] + isochrone.distance_modulus - mag_2_rgb = isochrone.mag_2[0: index_transition] + isochrone.distance_modulus - - mag_1_rgb = mag_1_rgb[::-1] - mag_2_rgb = mag_2_rgb[::-1] - - # Cut one way... - f_isochrone = scipy.interpolate.interp1d(mag_2_rgb, mag_1_rgb - mag_2_rgb, bounds_error=False, fill_value = 999.) - color_diff = np.fabs((g - r) - f_isochrone(r)) - cut_2 = (color_diff < np.sqrt(0.1**2 + r_err**2 + g_err**2)) - - # ...and now the other - f_isochrone = scipy.interpolate.interp1d(mag_1_rgb, mag_1_rgb - mag_2_rgb, bounds_error=False, fill_value = 999.) - color_diff = np.fabs((g - r) - f_isochrone(g)) - cut_1 = (color_diff < np.sqrt(0.1**2 + r_err**2 + g_err**2)) - - cut = np.logical_or(cut_1, cut_2) - - mag_bins = np.arange(17., 24.1, 0.1) - mag_centers = 0.5 * (mag_bins[1:] + mag_bins[0:-1]) - magerr = np.tile(0., len(mag_centers)) - for ii in range(0, len(mag_bins) - 1): - cut_mag_bin = (g > mag_bins[ii]) & (g < mag_bins[ii + 1]) - magerr[ii] = np.median(np.sqrt(0.1**2 + r_err[cut_mag_bin]**2 + g_err[cut_mag_bin]**2)) - - if return_all: - return cut, mag_centers[f_isochrone(mag_centers) < 100], (f_isochrone(mag_centers) + magerr)[f_isochrone(mag_centers) < 100], (f_isochrone(mag_centers) - magerr)[f_isochrone(mag_centers) < 100] - else: - return cut - -############################################################ - -def circle(x, y, r): - phi = np.linspace(0, 2 * np.pi, 1000) - return x + (r * np.cos(phi)), y + (r * np.sin(phi)) - -############################################################ - -try: - ra_select, dec_select = float(sys.argv[1]), float(sys.argv[2]) -except: - sys.exit('ERROR! Coordinates not given in correct format.') - -print('Search coordinates: (RA, Dec) = (%.2f, %.2f)'%(ra_select, dec_select)) - -# Now cut for a single pixel -pix_nside_select = ugali.utils.healpix.angToPix(nside, ra_select, dec_select) -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)]) - -############################################################ - -data_array = [] -for pix_nside in pix_nside_neighbors: - #infile = '%s/cat_hpx_%05i.fits'%(datadir, pix_nside) - infile = '%s/y3a2_ngmix_cm_%05i.fits'%(datadir, pix_nside) - #infile = '%s/*_%05i.fits'%(datadir, pix_nside) # TODO - get this to work for adaptable directories - if not os.path.exists(infile): - continue - reader = pyfits.open(infile) - data_array.append(reader[1].data) - reader.close() -print('Assembling data...') -data = np.concatenate(data_array) # TODO reduce this to just use needed columns so there is no excessive use of memory - -# De-redden magnitudes -#try: -# data = mlab.rec_append_fields(data, ['WAVG_MAG_PSF_DRED_G', 'WAVG_MAG_PSF_DRED_R'], [data[mag_g_dred_flag], data[mag_r_dred_flag]]) -#except: -# data = mlab.rec_append_fields(data, ['WAVG_MAG_PSF_DRED_G', 'WAVG_MAG_PSF_DRED_R'], [data[mag_g_flag] - data[extinction_g_flag], data[mag_r_flag] - data[extinction_r_flag]]) -#except: -# data = mlab.rec_append_fields(data, ['WAVG_MAG_PSF_DRED_G', 'WAVG_MAG_PSF_DRED_R'], [data[mag_g_flag], data[mag_r_flag]]) - -print('Found %i objects...'%(len(data))) - -############################################################ - -print('Applying cuts...') -#cut = (data[flags_g] < 4) & (data[flags_r] < 4) \ -# & (data['WAVG_MAG_PSF_DRED_G'] < 24.0) \ -# & ((data['WAVG_MAG_PSF_DRED_G'] - data['WAVG_MAG_PSF_DRED_R']) < 1.) \ -# & (np.fabs(data[spread_model_r_flag]) < 0.003 + data[spread_model_r_err_flag]) -# -#cut_gal = (data[flags_g] < 4) & (data[flags_r] < 4) \ -# & (data['WAVG_MAG_PSF_DRED_G'] < 24.0) \ -# & ((data['WAVG_MAG_PSF_DRED_G'] - data['WAVG_MAG_PSF_DRED_R']) < 1.) \ -# & (np.fabs(data[spread_model_r_flag]) > 0.003 + data[spread_model_r_err_flag]) - -#cut = (data[star_galaxy_classification] >= 0) & (data[star_galaxy_classification] <= 2) -#cut_gal = (data[star_galaxy_classification] > 2) -cut = (data[star_galaxy_classification] >= 0) & (data[star_galaxy_classification] <= 1) -cut_gal = (data[star_galaxy_classification] > 1) - -data_gal = data[cut_gal] -data = data[cut] - -print('%i star-like objects in ROI...'%(len(data))) -print('%i galaxy-like objects in ROI...'%(len(data_gal))) - -############################################################ - -if (cfg['fracdet'] is not None) and (cfg['fracdet'].lower().strip() != 'none') and (cfg['fracdet'] != ''): - print('Reading fracdet map %s ...'%(cfg['fracdet'])) - fracdet = healpy.read_map(cfg['fracdet']) -else: - print('No fracdet map specified ...') - fracdet = None - -############################################################ - -def searchByDistance(nside, data, distance_modulus, ra_select, dec_select, magnitude_threshold=24.5, plot=False, fracdet=None): - """ - Idea: - Send a data extension that goes to faint magnitudes, e.g., g < 24. - Use the whole region to identify hotspots using a slightly brighter - magnitude threshold, e.g., g < 23, so not susceptible to variations - in depth. Then compute the local field density using a small annulus - around each individual hotspot, e.g., radius 0.3 to 0.5 deg. - - plot = True enables diagnostic plotting for testings - fracdet corresponds to a fracdet map (numpy array, assumed to be EQUATORIAL and RING) - """ - - SCALE = 2.75 * (healpy.nside2resol(nside, arcmin=True) / 60.) # deg, scale for 2D histogram and various plotting - - print('Distance = %.1f kpc (m-M = %.1f)'%(ugali.utils.projector.distanceModulusToDistance(distance_modulus), distance_modulus)) - - dirname = '/home/s1/kadrlica/.ugali/isochrones/des/dotter2016/' - #dirname = '/Users/keithbechtol/Documents/DES/projects/mw_substructure/ugalidir/isochrones/des/dotter2016/' - iso = ugali.isochrone.factory('Dotter', hb_spread=0, dirname=dirname) - iso.age = 12. - iso.metallicity = 0.0001 - iso.distance_modulus = distance_modulus - - cut = cutIsochronePath(data[mag_g_dred_flag], data[mag_r_dred_flag], data[mag_err_g_flag], data[mag_err_r_flag], iso, radius=0.1) - data = data[cut] - cut_magnitude_threshold = (data[mag_g_dred_flag] < magnitude_threshold) - - print('%i objects left after isochrone cut...'%(len(data))) - - ### - - proj = ugali.utils.projector.Projector(ra_select, dec_select) - x, y = proj.sphereToImage(data['RA'][cut_magnitude_threshold], data['DEC'][cut_magnitude_threshold]) # Trimmed magnitude range for hotspot finding - x_full, y_full = proj.sphereToImage(data['RA'], data['DEC']) # In we want to use full magnitude range for significance evaluation - delta_x = 0.01 - area = delta_x**2 - smoothing = 2. / 60. # Was 3 arcmin - bins = np.arange(-1. * SCALE, SCALE + 1.e-10, delta_x) - centers = 0.5 * (bins[0: -1] + bins[1:]) - yy, xx = np.meshgrid(centers, centers) - - h = np.histogram2d(x, y, bins=[bins, bins])[0] - - if plot: - #pylab.figure('raw') - #pylab.clf() - #pylab.imshow(h.T, interpolation='nearest', extent=[-1. * SCALE, SCALE, -1. * SCALE, SCALE], origin='lower', cmap='binary') - #pylab.xlim([SCALE, -1. * SCALE]) - #pylab.ylim([-1. * SCALE, SCALE]) - #pylab.colorbar() - - ##import skymap - ##s = skymap.Skymap(projection='laea',llcrnrlon=340,llcrnrlat=-60,urcrnrlon=360,urcrnrlat=-50,lon_0 =355,lat_0=-55,celestial=False) - ##s.draw_hpxmap(fracdet) - #pix = skymap.healpix.ang2disc(355,-63,1) - #pix = skymap.healpix.ang2disc(4096,355,-63,1) - #s = skymap.Skymap() - #s.draw_hpxmap(m[pix],pix,4096) - #s.zoom_to_fit() - #s.zoom_to_fit(4096,m[pix],pix) - - reso = 0.25 - pylab.figure('gnom') - pylab.clf() - healpy.gnomview(fracdet, fig='gnom', rot=(ra_select, dec_select, 0.), reso=reso, xsize=(2. * SCALE * 60. / reso), - cmap='Greens', title='Fracdet') #binary - healpy.projscatter(data['RA'], data['DEC'], edgecolor='none', c='red', lonlat=True, s=2) - - h_g = scipy.ndimage.filters.gaussian_filter(h, smoothing / delta_x) - - #cut_goodcoverage = (data['NEPOCHS_G'][cut_magnitude_threshold] >= 2) & (data['NEPOCHS_R'][cut_magnitude_threshold] >= 2) - # expect NEPOCHS to be good in DES data - - delta_x_coverage = 0.1 - area_coverage = (delta_x_coverage)**2 - bins_coverage = np.arange(-5., 5. + 1.e-10, delta_x_coverage) - h_coverage = np.histogram2d(x, y, bins=[bins_coverage, bins_coverage])[0] - #h_goodcoverage = np.histogram2d(x[cut_goodcoverage], y[cut_goodcoverage], bins=[bins_coverage, bins_coverage])[0] - h_goodcoverage = np.histogram2d(x, y, bins=[bins_coverage, bins_coverage])[0] - - n_goodcoverage = h_coverage[h_goodcoverage > 0].flatten() - - #characteristic_density = np.mean(n_goodcoverage) / area_coverage # per square degree - characteristic_density = np.median(n_goodcoverage) / area_coverage # per square degree - print('Characteristic density = %.1f deg^-2'%(characteristic_density)) - - # Use pixels with fracdet ~1.0 to estimate the characteristic density - if fracdet is not None: - fracdet_zero = np.tile(0., len(fracdet)) - cut = (fracdet != healpy.UNSEEN) - fracdet_zero[cut] = fracdet[cut] - - nside_fracdet = healpy.npix2nside(len(fracdet)) - - subpix_region_array = [] - for pix in np.unique(ugali.utils.healpix.angToPix(nside, data['RA'], data['DEC'])): - subpix_region_array.append(ugali.utils.healpix.subpixel(pix, nside, nside_fracdet)) - subpix_region_array = np.concatenate(subpix_region_array) - - # Compute mean fracdet in the region so that this is available as a correction factor - cut = (fracdet[subpix_region_array] != healpy.UNSEEN) - mean_fracdet = np.mean(fracdet[subpix_region_array[cut]]) - - subpix_region_array = subpix_region_array[fracdet[subpix_region_array] > 0.99] - subpix = ugali.utils.healpix.angToPix(nside_fracdet, - data['RA'][cut_magnitude_threshold], - data['DEC'][cut_magnitude_threshold]) # Remember to apply mag threshold to objects - characteristic_density_fracdet = float(np.sum(np.in1d(subpix, subpix_region_array))) \ - / (healpy.nside2pixarea(nside_fracdet, degrees=True) * len(subpix_region_array)) # deg^-2 - print('Characteristic density fracdet = %.1f deg^-2'%(characteristic_density_fracdet)) - - # Correct the characteristic density by the mean fracdet value - characteristic_density_raw = 1. * characteristic_density - characteristic_density /= mean_fracdet - print('Characteristic density (fracdet corrected) = %.1f deg^-2'%(characteristic_density)) - - if plot: - pylab.figure('poisson') - pylab.clf() - n_max = np.max(h_coverage) - pylab.hist(n_goodcoverage, bins=np.arange(n_max) - 0.5, color='blue', histtype='step', lw=2, normed=True) - pylab.scatter(np.arange(n_max), scipy.stats.poisson.pmf(np.arange(n_max), mu=np.median(n_goodcoverage)), c='red', edgecolor='none', zorder=10) - #pylab.plot(np.arange(n_max), scipy.stats.poisson.pmf(np.arange(n_max), mu=np.median(n_goodcoverage)), c='red', lw=2, zorder=10) - pylab.axvline(characteristic_density * area_coverage, c='black', ls='--') - if fracdet is not None: - pylab.axvline(characteristic_density_raw * area_coverage, c='orange', ls='--') - pylab.axvline(characteristic_density_fracdet * area_coverage, c='green', ls='--') - pylab.xlabel('Counts per %.3f deg^-2 pixel'%(area_coverage)) - pylab.ylabel('PDF') - - if plot: - vmax = min(3. * characteristic_density * area, np.max(h_g)) - - pylab.figure('smooth') - pylab.clf() - pylab.imshow(h_g.T, - interpolation='nearest', extent=[-1. * SCALE, SCALE, -1. * SCALE, SCALE], origin='lower', cmap='gist_heat', vmax=vmax) - pylab.colorbar() - pylab.xlim([SCALE, -1. * SCALE]) - pylab.ylim([-1. * SCALE, SCALE]) - pylab.xlabel(r'$\Delta$ RA (deg)') - pylab.ylabel(r'$\Delta$ Dec (deg)') - pylab.title('(RA, Dec, mu) = (%.2f, %.2f, %.2f)'%(ra_select, dec_select, distance_modulus)) - - factor_array = np.arange(1., 5., 0.05) - rara, decdec = proj.imageToSphere(xx.flatten(), yy.flatten()) - cutcut = (ugali.utils.healpix.angToPix(nside, rara, decdec) == pix_nside_select).reshape(xx.shape) - threshold_density = 5 * characteristic_density * area - for factor in factor_array: - h_region, n_region = scipy.ndimage.measurements.label((h_g * cutcut) > (area * characteristic_density * factor)) - #print 'factor', factor, n_region, n_region < 10 - if n_region < 10: - threshold_density = area * characteristic_density * factor - break - - h_region, n_region = scipy.ndimage.measurements.label((h_g * cutcut) > threshold_density) - h_region = np.ma.array(h_region, mask=(h_region < 1)) - if plot: - pylab.figure('regions') - pylab.clf() - pylab.imshow(h_region.T, - interpolation='nearest', extent=[-1. * SCALE, SCALE, -1. * SCALE, SCALE], origin='lower') - pylab.colorbar() - pylab.xlim([SCALE, -1. * SCALE]) - pylab.ylim([-1. * SCALE, SCALE]) - - ra_peak_array = [] - dec_peak_array = [] - r_peak_array = [] - sig_peak_array = [] - distance_modulus_array = [] - - pylab.figure('sig') - pylab.clf() - for index in range(1, n_region + 1): - index_peak = np.argmax(h_g * (h_region == index)) - x_peak, y_peak = xx.flatten()[index_peak], yy.flatten()[index_peak] - #print index, np.max(h_g * (h_region == index)) - if plot: - pylab.figure('regions') - #pylab.scatter(x_peak, y_peak, marker='x', c='white') - pylab.scatter(x_peak, y_peak, marker='o', c='none', edgecolor='black', s=50) - - - #angsep_peak = np.sqrt((x_full - x_peak)**2 + (y_full - y_peak)**2) # Use full magnitude range, NOT TESTED!!! - angsep_peak = np.sqrt((x - x_peak)**2 + (y - y_peak)**2) # Impose magnitude threshold - - # Compute the local characteristic density - # If fracdet map is available, use that information to either compute local density, - # or in regions of spotty coverage, use the typical density of the region - if fracdet is not None: - ra_peak, dec_peak = proj.imageToSphere(x_peak, y_peak) - subpix_all = ugali.utils.healpix.angToDisc(nside_fracdet, ra_peak, dec_peak, 0.5) - subpix_inner = ugali.utils.healpix.angToDisc(nside_fracdet, ra_peak, dec_peak, 0.3) - subpix_annulus = subpix_all[~np.in1d(subpix_all, subpix_inner)] - mean_fracdet = np.mean(fracdet_zero[subpix_annulus]) - print('mean_fracdet', mean_fracdet) - if mean_fracdet < 0.5: - characteristic_density_local = characteristic_density - print('characteristic_density_local baseline', characteristic_density_local) - else: - # Check pixels in annulus with complete coverage - subpix_annulus_region = np.intersect1d(subpix_region_array, subpix_annulus) - print(float(len(subpix_annulus_region)) / len(subpix_annulus)) - if (float(len(subpix_annulus_region)) / len(subpix_annulus)) < 0.25: - characteristic_density_local = characteristic_density - print('characteristic_density_local spotty', characteristic_density_local) - else: - characteristic_density_local = float(np.sum(np.in1d(subpix, subpix_annulus_region))) \ - / (healpy.nside2pixarea(nside_fracdet, degrees=True) * len(subpix_annulus_region)) # deg^-2 - print('characteristic_density_local cleaned up', characteristic_density_local) - else: - # Compute the local characteristic density - area_field = np.pi * (0.5**2 - 0.3**2) - n_field = np.sum((angsep_peak > 0.3) & (angsep_peak < 0.5)) - characteristic_density_local = n_field / area_field - - # If not good azimuthal coverage, revert - cut_annulus = (angsep_peak > 0.3) & (angsep_peak < 0.5) - #phi = np.degrees(np.arctan2(y_full[cut_annulus] - y_peak, x_full[cut_annulus] - x_peak)) # Use full magnitude range, NOT TESTED!!! - phi = np.degrees(np.arctan2(y[cut_annulus] - y_peak, x[cut_annulus] - x_peak)) # Impose magnitude threshold - h = np.histogram(phi, bins=np.linspace(-180., 180., 13))[0] - if np.sum(h > 0) < 10 or np.sum(h > 0.5 * np.median(h)) < 10: - #angsep_peak = np.sqrt((x - x_peak)**2 + (y - y_peak)**2) - characteristic_density_local = characteristic_density - - print('Characteristic density local = %.1f deg^-2'%(characteristic_density_local)) - - size_array = np.arange(0.01, 0.3, 0.01) - sig_array = np.tile(0., len(size_array)) - for ii in range(0, len(size_array)): - n_peak = np.sum(angsep_peak < size_array[ii]) - n_model = characteristic_density_local * (np.pi * size_array[ii]**2) - sig_array[ii] = scipy.stats.norm.isf(scipy.stats.poisson.sf(n_peak, n_model)) - if sig_array[ii] > 25: - sig_array[ii] = 25. # Set a maximum significance value - - if plot: - pylab.figure('sig') - pylab.plot(size_array, sig_array) - pylab.xlabel('Radius used to compute significance (deg)') - pylab.ylabel('Detection significance') - - ra_peak, dec_peak = proj.imageToSphere(x_peak, y_peak) - r_peak = size_array[np.argmax(sig_array)] - if np.max(sig_array) == 25.: - r_peak = 0.5 - - #print 'Candidate:', x_peak, y_peak, r_peak, np.max(sig_array), ra_peak, dec_peak - print('Candidate: %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f'%(x_peak, y_peak, r_peak, np.max(sig_array), ra_peak, dec_peak)) - if np.max(sig_array) > 5.: - if plot: - x_circle, y_circle = circle(x_peak, y_peak, r_peak) - pylab.figure('smooth') - pylab.plot(x_circle, y_circle, c='gray') - pylab.text(x_peak - r_peak, y_peak + r_peak, r'%.2f $\sigma$'%(np.max(sig_array)), color='gray') - - ra_peak_array.append(ra_peak) - dec_peak_array.append(dec_peak) - r_peak_array.append(r_peak) - sig_peak_array.append(np.max(sig_array)) - distance_modulus_array.append(distance_modulus) - - if plot: - input('Plots are ready...') - - return ra_peak_array, dec_peak_array, r_peak_array, sig_peak_array, distance_modulus_array - -############################################################ - -def diagnostic(data, data_gal, ra_peak, dec_peak, r_peak, sig_peak, distance_modulus, age=12., metallicity=0.0001): - - # Dotter isochrones - dirname = '/home/s1/kadrlica/.ugali/isochrones/des/dotter2016/' - #dirname = '/Users/keithbechtol/Documents/DES/projects/mw_substructure/ugalidir/isochrones/des/dotter2016/' - iso = ugali.isochrone.factory('Dotter', hb_spread=0, dirname=dirname) - iso.age = age - iso.metallicity = metallicity - iso.distance_modulus = distance_modulus - - ## Padova isochrones - #dirname_alt = '/home/s1/kadrlica/.ugali/isochrones/des/bressan2012/' #padova/' - #iso_alt = ugali.isochrone.factory('Padova', hb_spread=0, dirname=dirname_alt) - #iso_alt.age = age - #iso_alt.metallicity = metallicity - #iso_alt.distance_modulus = distance_modulus - - cut_iso, g_iso, gr_iso_min, gr_iso_max = cutIsochronePath(data[mag_g_dred_flag], - data[mag_r_dred_flag], - data[mag_err_g_flag], - data[mag_err_r_flag], - iso, - radius=0.1, - return_all=True) - - #cut_iso_gal = cutIsochronePath(data_gal['MAG_PSF_G'], # TODO: should these be WAVG? and if so, should these then also be dereddened flags? - # data_gal['MAG_PSF_R'], - # data_gal['MAGERR_PSF_G'], - # data_gal['MAGERR_PSF_R'], - # iso, - # radius=0.1, - # return_all=False) - cut_iso_gal = cutIsochronePath(data_gal[mag_g_dred_flag], - data_gal[mag_r_dred_flag], - data_gal[mag_err_g_flag], - data_gal[mag_err_r_flag], - iso, - radius=0.1, - return_all=False) - - proj = ugali.utils.projector.Projector(ra_peak, dec_peak) - x, y = proj.sphereToImage(data['RA'][cut_iso], data['DEC'][cut_iso]) - x_gal, y_gal = proj.sphereToImage(data_gal['RA'][cut_iso_gal], data_gal['DEC'][cut_iso_gal]) - -########################################################### - - angsep = ugali.utils.projector.angsep(ra_peak, dec_peak, data['RA'], data['DEC']) - cut_inner = (angsep < r_peak) - cut_annulus = (angsep > 0.5) & (angsep < 1.) - - angsep_gal = ugali.utils.projector.angsep(ra_peak, dec_peak, data_gal['RA'], data_gal['DEC']) - cut_inner_gal = (angsep_gal < r_peak) - cut_annulus_gal = (angsep_gal > 0.5) & (angsep_gal < 1.) - -# ########## -# -# # Check for possible associations -# glon_peak, glat_peak = ugali.utils.projector.celToGal(ra_peak, dec_peak) -# catalog_array = ['McConnachie12', '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 = '; %s at %.3f deg'%(match[0]['name'], sep) -# else: -# association_string = '; no association within 0.5 deg' - -############################################################ - -distance_modulus_search_array = np.arange(16., 24., 0.5) - -ra_peak_array = [] -dec_peak_array = [] -r_peak_array = [] -sig_peak_array = [] -distance_modulus_array = [] -for distance_modulus in distance_modulus_search_array: - ra_peak, dec_peak, r_peak, sig_peak, distance_modulus = searchByDistance(nside, data, distance_modulus, ra_select, dec_select, fracdet=fracdet) - ra_peak_array.append(ra_peak) - dec_peak_array.append(dec_peak) - r_peak_array.append(r_peak) - sig_peak_array.append(sig_peak) - distance_modulus_array.append(distance_modulus) - -ra_peak_array = np.concatenate(ra_peak_array) -dec_peak_array = np.concatenate(dec_peak_array) -r_peak_array = np.concatenate(r_peak_array) -sig_peak_array = np.concatenate(sig_peak_array) -distance_modulus_array = np.concatenate(distance_modulus_array) - -# Sort peaks according to significance -index_sort = np.argsort(sig_peak_array)[::-1] -ra_peak_array = ra_peak_array[index_sort] -dec_peak_array = dec_peak_array[index_sort] -r_peak_array = r_peak_array[index_sort] -sig_peak_array = sig_peak_array[index_sort] -distance_modulus_array = distance_modulus_array[index_sort] - -# Collect overlapping peaks -for ii in range(0, len(sig_peak_array)): - if sig_peak_array[ii] < 0: - continue - angsep = ugali.utils.projector.angsep(ra_peak_array[ii], dec_peak_array[ii], ra_peak_array, dec_peak_array) - sig_peak_array[(angsep < r_peak_array[ii]) & (np.arange(len(sig_peak_array)) > ii)] = -1. - -# Prune the list of peaks -ra_peak_array = ra_peak_array[sig_peak_array > 0.] -dec_peak_array = dec_peak_array[sig_peak_array > 0.] -r_peak_array = r_peak_array[sig_peak_array > 0.] -distance_modulus_array = distance_modulus_array[sig_peak_array > 0.] -sig_peak_array = sig_peak_array[sig_peak_array > 0.] # Update the sig_peak_array last! - -for ii in range(0, len(sig_peak_array)): - print('%.2f sigma; (RA, Dec, d) = (%.2f, %.2f); r = %.2f deg; d = %.1f, mu = %.2f mag)'%(sig_peak_array[ii], - ra_peak_array[ii], - dec_peak_array[ii], - r_peak_array[ii], - ugali.utils.projector.distanceModulusToDistance(distance_modulus_array[ii]), - distance_modulus_array[ii])) - - if (sig_peak_array[ii] > 5.5) & (r_peak_array[ii] < 0.28): - diagnostic(data, data_gal, ra_peak_array[ii], dec_peak_array[ii], r_peak_array[ii], sig_peak_array[ii], distance_modulus_array[ii]) - - -outfile = '%s/results_nside_%s_%i.csv'%(results_dir, nside, pix_nside_select) -writer = open(outfile, 'w') -#writer.write('sig, ra, dec, distance_modulus, r\n') -for ii in range(0, len(sig_peak_array)): - writer.write('%10.2f, %10.2f, %10.2f, %10.2f, %10.2f\n'%(sig_peak_array[ii], - ra_peak_array[ii], - dec_peak_array[ii], - distance_modulus_array[ii], - r_peak_array[ii])) -writer.close()