In [ ]:
import os
import pickle
import sys
globber_path = '/Users/adrian/projects/globber/'
if globber_path not in sys.path:
sys.path.append(globber_path)
# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.table as table
import astropy.units as u
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as pl
pl.style.use('apw-notebook')
%matplotlib inline
from astropy.io.fits import getdata
from scipy.misc import logsumexp
from scipy.ndimage import gaussian_filter
from scipy import interpolate
from scipy.stats import scoreatpercentile
from astroML.utils import log_multivariate_gaussian
import h5py
from globber.ngc5897 import cluster_c, r_c, r_t
In [ ]:
XCov_filename = "/Users/adrian/projects/globber/data/ngc5897/XCov_lg.h5"
In [ ]:
with h5py.File(XCov_filename, "r") as f:
allX = f['search']['X'][:]
pre_filter_ix = (allX[:,0] > 18.) & (allX[:,0] < 21.) & (allX[:,1] > 0.) & (allX[:,1] < 1.)
allX = allX[pre_filter_ix]
ra = f['search']['ra'][:][pre_filter_ix]
dec = f['search']['dec'][:][pre_filter_ix]
cluX = f['cluster']['X'][:]
all_c = coord.SkyCoord(ra=ra*u.degree, dec=dec*u.degree)
In [ ]:
pl.figure(figsize=(8,6))
pl.plot(ra, dec, ls='none', marker=',', alpha=0.05)
pl.gca().set_aspect('equal')
pl.xlim(pl.xlim()[::-1])
In [ ]:
cluster_ix = all_c.separation(cluster_c) < 6*u.arcmin
In [ ]:
def search_field(ra, dec):
ix1 = (ra > 215) & (ra < 240) & (dec < -17) & (dec > -30)
return ix1
search_ix = search_field(ra, dec)
In [ ]:
# def control_field(ra, dec):
# ix1 = (ra > 228) & (ra < 232) & (dec < -24) & (dec > -26)
# ix2 = (ra > 228) & (ra < 232) & (dec < -17) & (dec > -19)
# return ix1 | ix2
# control_ix = control_field(ra, dec)
control_ix = search_ix & np.logical_not(cluster_ix)
In [ ]:
pl.figure(figsize=(8,6))
pl.plot(ra, dec, ls='none', marker=',', alpha=0.1)
pl.plot(ra[search_ix], dec[search_ix], ls='none', marker=',', alpha=0.1, color='g')
pl.plot(ra[control_ix], dec[control_ix], ls='none', marker=',', alpha=0.1, color='r')
pl.gca().set_aspect('equal')
pl.xlim(pl.xlim()[::-1])
# pl.xlim(235,220)
# pl.ylim(-26,-16)
In [ ]:
searchX = allX[search_ix]
TODO: a script that does this for steps in dm of 0.05 from 0 to 1 for normed and un-normed
In [ ]:
dm_offset = 0.
In [ ]:
search_data = dict()
cluster_data = dict()
control_data = dict()
color_names = ['g-r','g-i','g-z']
names = ['i'] + color_names
for i,name in zip([0,1,2,3], names):
search_data[name] = searchX[:,i]
cluster_data[name] = cluX[:,i]
control_data[name] = allX[control_ix,i]
search_data = table.Table(search_data).as_array()
cluster_data = table.Table(cluster_data).as_array()
control_data = table.Table(control_data).as_array()
cluster_data['i'] += dm_offset
In [ ]:
c_step = 0.03
m_step = 0.06
c_bins = np.arange(0,0.7+c_step,c_step)
m_bins = np.arange(18,21+m_step,m_step)
In [ ]:
def bin_cmd(color, mag):
H,_,_ = np.histogram2d(color, mag, bins=(c_bins, m_bins))
m_mesh,c_mesh = np.meshgrid((m_bins[1:]+m_bins[:-1])/2, (c_bins[1:]+c_bins[:-1])/2)
return c_mesh, m_mesh, H
In [ ]:
cluster_Hs = dict()
control_Hs = dict()
control_spline_Hs = dict()
for cname in color_names:
xmesh,ymesh,cluster_H = bin_cmd(cluster_data[cname], cluster_data['i'])
xmesh,ymesh,control_H = bin_cmd(control_data[cname], control_data['i'])
cluster_Hs[cname] = cluster_H / cluster_H.sum()
control_Hs[cname] = control_H / control_H.sum()
# use a bivariate spline to smooth the control CMD
spl = interpolate.SmoothBivariateSpline(xmesh.ravel(), ymesh.ravel(),
control_H.ravel()/control_H.sum(), kx=5, ky=5)
spl_control_H = spl.ev(xmesh.ravel(), ymesh.ravel())
spl_control_H = spl_control_H.reshape(xmesh.shape)
control_spline_Hs[cname] = spl_control_H
# ---------------------------------------------
# diagnostic plot
fig,axes = pl.subplots(1,3,figsize=(10,6),sharex=True,sharey=True)
ax = axes[0]
ax.pcolormesh(xmesh, ymesh, cluster_H, cmap='Blues')
ax.set_xlim(xmesh.min(), xmesh.max())
ax.set_ylim(ymesh.max(), ymesh.min())
ax.set_xlabel('${}$'.format(cname))
ax.set_ylabel('$i$')
ax.set_title("Cluster stars")
ax = axes[1]
ax.pcolormesh(xmesh, ymesh, control_H, cmap='Blues')
ax.set_xlabel('${}$'.format(cname))
ax.set_title("Control field")
ax = axes[2]
ax.pcolormesh(xmesh, ymesh, spl_control_H, cmap='Blues')
ax.set_xlabel('${}$'.format(cname))
ax.set_title("Spline smoothed")
In [ ]:
fig,axes = pl.subplots(1,len(color_names),figsize=(12,4),sharex=True,sharey=True)
for i,cname in enumerate(color_names):
ax = axes[i]
_,bins,_ = ax.hist(cluster_Hs[cname].ravel(), bins=np.logspace(-5,-1,32), alpha=0.4)
ax.hist(control_spline_Hs[cname].ravel(), bins=bins, alpha=0.4)
ax.set_xscale('log')
ax.set_yscale('log')
In [ ]:
fig,axes = pl.subplots(1,len(color_names),figsize=(12,4),sharex=True,sharey=True)
for i,cname in enumerate(color_names):
ax = axes[i]
div = cluster_Hs[cname] / control_spline_Hs[cname]
div[div < 0.] = 0.
ax.hist(div.ravel(), bins=np.logspace(-1.5,2.5,32), alpha=0.4)
ax.set_xscale('log')
ax.set_yscale('log')
In [ ]:
# thresh = 20
# null_thresh = 4E-4
threshs = {
'g-r': 10,
'g-i': 15,
'g-z': 25,
}
null_threshs = {
'g-r': 25E-4,
'g-i': 20E-4,
'g-z': 2E-3,
}
matched_filters = dict()
for cname in color_names:
null_thresh = null_threshs[cname]
thresh = threshs[cname]
div = cluster_Hs[cname] / control_spline_Hs[cname]
div[div < 0.] = 0.
div[div > thresh] = thresh
div[(cluster_Hs[cname] < null_thresh) & (control_spline_Hs[cname] < null_thresh)] = 0.
matched_filters[cname] = div.copy()
fig,axes = pl.subplots(1,3,figsize=(10,6),sharex=True,sharey=True)
ax = axes[0]
ax.pcolormesh(xmesh, ymesh, cluster_Hs[cname], cmap='Blues')
ax.set_xlim(xmesh.min(), xmesh.max())
ax.set_ylim(ymesh.max(), ymesh.min())
ax.set_xlabel('${}$'.format(cname))
ax.set_ylabel('$i$')
ax.set_title("Cluster stars")
ax = axes[1]
ax.pcolormesh(xmesh, ymesh, spl_control_H, cmap='Blues')
ax.set_xlabel('${}$'.format(cname))
ax.set_title("Control field")
ax = axes[2]
ax.pcolormesh(xmesh, ymesh, div, cmap='Blues')
ax.set_xlabel('${}$'.format(cname))
ax.set_title("Spline smoothed")
In [ ]:
cmd_smooth = 0.02 # mag
smooth_matched_filters = dict()
for cname in color_names:
smooth_matched_filters[cname] = gaussian_filter(matched_filters[cname],
sigma=[cmd_smooth/c_step,cmd_smooth/m_step])
fig,axes = pl.subplots(1,len(color_names),figsize=(10,6),sharex=True,sharey=True)
for i,cname in enumerate(color_names):
ax = axes[i]
ax.pcolormesh(xmesh, ymesh, smooth_matched_filters[cname], cmap='Blues')
ax.set_xlabel('${}$'.format(cname))
ax.set_xlim(xmesh.min(), xmesh.max())
ax.set_ylim(ymesh.max(), ymesh.min())
In [ ]:
print(search_ix.sum())
In [ ]:
n_search = search_ix.sum()
weights = dict()
for cname in color_names:
# for cname in ['g-z']:
_x,_y = digitize2d(search_data[cname], search_data['i'], c_bins, m_bins)
weights[cname] = np.zeros(n_search)
for i in range(n_search):
try:
weights[cname][i] = smooth_matched_filters[cname][_x[i]-1,_y[i]-1]
except IndexError:
weights[cname][i] = 0. # falls outside bin??
# weights[cname] = np.array([smooth_matched_filters[cname][_x[i]-1,_y[i]-1] for i in range(n_search)])
In [ ]:
# n_search = search_ix.sum()
# weights = dict()
# for cname in color_names:
# weights[cname] = np.zeros(n_search)
# for i in range(n_search):
# for j,cname in enumerate(color_names):
# derp,_,_ = np.histogram2d(search_data[cname][i:i+1], search_data['i'][i:i+1], bins=(c_bins, m_bins))
# weights[cname][i] = (derp * smooth_matched_filters[cname]).sum()
# if (i % 10000) == 0:
# print(i)
In [ ]:
sky_binsize = (6*u.arcmin).to(u.degree).value
sky_smooth = (6*u.arcmin).to(u.degree).value / sky_binsize
search_ra = ra[search_ix]
search_dec = dec[search_ix]
ra_bins = np.arange(search_ra.min(), search_ra.max()+sky_binsize, sky_binsize)
dec_bins = np.arange(search_dec.min(), search_dec.max()+sky_binsize, sky_binsize)
In [ ]:
search_H_sky = None
unw_search_H_sky = None
for cname in weights.keys():
_H_sky,ra_edges,dec_edges = np.histogram2d(search_ra, search_dec,
bins=(ra_bins, dec_bins),
weights=weights[cname])
_unw_H_sky,ra_edges,dec_edges = np.histogram2d(search_ra, search_dec,
bins=(ra_bins, dec_bins))
if search_H_sky is None:
search_H_sky = _H_sky.T
unw_search_H_sky = _unw_H_sky.T
else:
# search_H_sky += _H_sky.T
# unw_search_H_sky += _unw_H_sky.T
search_H_sky *= _H_sky.T
unw_search_H_sky *= _unw_H_sky.T
norm_search_H_sky = search_H_sky / unw_search_H_sky
ra_mesh,dec_mesh = np.meshgrid((ra_edges[1:]+ra_edges[:-1])/2, (dec_edges[1:]+dec_edges[:-1])/2)
In [ ]:
# save_path = "/Users/adrian/projects/globber/data/ngc5897/density_maps"
# if not os.path.exists(save_path):
# os.mkdir(save_path)
# np.save(os.path.join(save_path, 'ddm{:.2f}.npy'.format(dm_offset)), search_H_sky)
# np.save(os.path.join(save_path, 'ddm{:.2f}_norm.npy'.format(dm_offset)), norm_search_H_sky)
In [ ]:
# H_operation = lambda x: np.log(x)
H_operation = lambda x: np.sqrt(x)
# H_operation = lambda x: x
# H_operation = lambda x: x**2
In [ ]:
tmp = H_operation(norm_search_H_sky.ravel())
bins = np.linspace(*scoreatpercentile(tmp, [1,99]), num=32)
pl.hist(tmp, bins=bins);
vmin,vmax = scoreatpercentile(tmp, [15,85])
pl.axvline(vmin, color='r')
pl.axvline(vmax, color='r')
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(15,6),sharex=True,sharey=True)
ax = axes[0]
ax.pcolormesh(ra_mesh, dec_mesh, H_operation(norm_search_H_sky),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_xlim(ra_mesh.max(), ra_mesh.min())
ax.set_ylim(dec_mesh.min(), dec_mesh.max())
ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_aspect('equal')
ax = axes[1]
ax.pcolormesh(ra_mesh, dec_mesh, gaussian_filter(H_operation(norm_search_H_sky), sky_smooth),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_xlim(ra_mesh.max(), ra_mesh.min())
ax.set_ylim(dec_mesh.min(), dec_mesh.max())
ax.set_xlabel('RA [deg]')
ax.set_aspect('equal')
In [ ]:
zoom_buffer = 3
fig,axes = pl.subplots(1,2,figsize=(15,6),sharex=True,sharey=True)
ax = axes[0]
ax.pcolormesh(ra_mesh, dec_mesh, H_operation(norm_search_H_sky),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_xlim(cluster_c.ra.degree+zoom_buffer, cluster_c.ra.degree-zoom_buffer)
ax.set_ylim(cluster_c.dec.degree-zoom_buffer, cluster_c.dec.degree+zoom_buffer)
ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_aspect('equal')
ax = axes[1]
ax.pcolormesh(ra_mesh, dec_mesh, gaussian_filter(H_operation(norm_search_H_sky), sky_smooth),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_xlim(cluster_c.ra.degree+zoom_buffer, cluster_c.ra.degree-zoom_buffer)
ax.set_ylim(cluster_c.dec.degree-zoom_buffer, cluster_c.dec.degree+zoom_buffer)
ax.set_xlabel('RA [deg]')
ax.set_aspect('equal')
In [ ]:
tmp = np.log10(search_H_sky.ravel())
fig,ax = pl.subplots(1,1,figsize=(6,6),sharex=True,sharey=True)
# ax.contour(ra_mesh, dec_mesh, search_H_sky, levels=10**np.linspace(nvmin,tmp.max(),12), colors='k') # cmap='magma_r',
ax.contour(ra_mesh, dec_mesh, search_H_sky, levels=np.logspace(np.log10(vmin),tmp.max(),10), cmap='magma')
# ax.set_xlim(ra_mesh.max(), ra_mesh.min())
# ax.set_ylim(dec_mesh.min(), dec_mesh.max())
ax.set_xlim(cluster_c.ra.degree + zoom_buffer, cluster_c.ra.degree - zoom_buffer)
ax.set_ylim(cluster_c.dec.degree - zoom_buffer, cluster_c.dec.degree + zoom_buffer)
ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_aspect('equal')
In [ ]:
# # H_operation = lambda x: np.log(x)
# # H_operation = lambda x: np.sqrt(x)
# H_operation = lambda x: x
# # H_operation = lambda x: x**2
In [ ]:
# zoom_buffer = 5
# vmin = None
# vmax = None
# for j,dm_offset in enumerate([-0.4,-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.4]):
# the_H = np.load(os.path.join(save_path, 'ddm{:.2f}.npy'.format(dm_offset)))
# if vmin is None:
# tmp = H_operation(the_H.ravel())
# vmin,vmax = scoreatpercentile(tmp, [15,85])
# fig,ax = pl.subplots(1,1,figsize=(6,6))
# ax.pcolormesh(ra_mesh, dec_mesh, H_operation(the_H),
# cmap='Greys', vmin=vmin, vmax=vmax)
# pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
# edgecolor='r', facecolor='none')
# ax.add_patch(pa)
# pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
# edgecolor='g', facecolor='none')
# ax.add_patch(pa)
# ax.set_xlim(cluster_c.ra.degree+zoom_buffer, cluster_c.ra.degree-zoom_buffer)
# ax.set_ylim(cluster_c.dec.degree-zoom_buffer, cluster_c.dec.degree+zoom_buffer)
# ax.set_xlabel('RA [deg]')
# ax.set_ylabel('Dec [deg]')
# ax.set_aspect('equal')
# fig.savefig("/Users/adrian/Downloads/{}.png".format(j), dpi=300)