From Brani's PS1 query, dereddened with scripts/deredden-ps1.py
In [ ]:
import os
import pickle
import sys
# 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
import healpy as hp
from scipy.misc import logsumexp
from scipy.ndimage import gaussian_filter
from astroML.utils import log_multivariate_gaussian
import gary.coordinates as gc
import gary.dynamics as gd
import gary.potential as gp
from gary.observation import distance_modulus
from ophiuchus import galactocentric_frame, vcirc, vlsr
import ophiuchus.potential as op
In [ ]:
catalog_data_path = "/Users/adrian/projects/globber/data/gc_catalogs/"
data_path = "/Users/adrian/projects/globber/data/ngc5897/"
plot_path = "/Users/adrian/projects/globber/plots/ngc5897"
if not os.path.exists(plot_path):
os.mkdir(plot_path)
In [ ]:
cluster_name = "NGC 5897"
# Global configuration stuff
DM = 15.62 # Brani's fit
cluster_c = coord.SkyCoord(ra=229.352*u.degree,
dec=-21.01*u.degree) # cluster sky position
in_cluster_pad = 0.075*u.degree # include cluster stars within (small for purity)
out_cluster_pad = 0.12*u.degree # exclude cluster stars within (large to be safe)
g_lim = (17.5, 20.5) # range of g-band magnitudes to use
color_lim = {
'g-r': (0,0.5),
'g-i': (0,0.75),
'g-z': (0,0.75),
}
In [ ]:
pm_gc_main = np.genfromtxt(os.path.join(catalog_data_path,"gl_2012_J2000.cat1.txt"), dtype=None,
skip_header=2,
usecols=[0,2,3,6,7,8,9,10,11,12,13],
names=['ngc_num','ra','dec','dist','dist_err','mu_ra','mu_ra_err',
'mu_dec','mu_dec_err', 'vr', 'vr_err'])
pm_gc_main = table.Table(pm_gc_main)
In [ ]:
go = ascii.read(os.path.join(catalog_data_path,"go97_table1.txt"))
In [ ]:
all_gc = pm_gc_main
all_gc['name'] = np.array(["NGC {}".format(x) for x in all_gc['ngc_num']])
all_gc = table.join(all_gc, go, keys='name')
cluster = all_gc[all_gc['name'] == cluster_name]
# HACK: use Branimir's distance modulus and make uncertainty better
from gary.observation import distance as distance_from_dm
cluster['dist'] = distance_from_dm(15.62).to(u.kpc).value
cluster['dist_err'] = 0.05*cluster['dist']
cluster
In [ ]:
ps1 = np.load(os.path.join(data_path, "PS1_stars_pv3_dered_sm.npy"))
ps1 = ps1[np.isfinite(ps1['dered_z'])]
ps1_c = coord.ICRS(ra=ps1['ra']*u.degree, dec=ps1['dec']*u.degree)
In [ ]:
ps1_ra_bounds = (ps1_c.ra.min(), ps1_c.ra.max())
ps1_dec_bounds = (ps1_c.dec.min(), ps1_c.dec.max())
ps1_ra_bounds, ps1_dec_bounds
In [ ]:
nsamples = 128
In [ ]:
print(cluster.colnames)
In [ ]:
cluster_c = coord.ICRS(ra=cluster['ra']*u.degree,
dec=cluster['dec']*u.degree,
distance=cluster['dist']*u.kpc)
In [ ]:
xyz = cluster_c.transform_to(galactocentric_frame).cartesian.xyz
vxyz = gc.vhel_to_gal(cluster_c, pm=(cluster['mu_ra']*u.mas/u.yr,
cluster['mu_dec']*u.mas/u.yr),
rv=cluster['vr']*u.km/u.s,
galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
In [ ]:
np.random.seed(42)
_distances = np.random.normal(cluster['dist'], cluster['dist_err'], size=nsamples)
cluster_samples_c = coord.ICRS(ra=(np.zeros(nsamples) + cluster['ra'])*u.degree,
dec=(np.zeros(nsamples) + cluster['dec'])*u.degree,
distance=_distances*u.kpc)
_mu_ras = np.random.normal(cluster['mu_ra'], cluster['mu_ra_err'], size=nsamples)
_mu_decs = np.random.normal(cluster['mu_dec'], cluster['mu_dec_err'], size=nsamples)
_vrs = np.random.normal(cluster['vr'], cluster['vr_err'], size=nsamples)
# ---
samples_xyz = cluster_samples_c.transform_to(galactocentric_frame).cartesian.xyz
samples_vxyz = gc.vhel_to_gal(cluster_samples_c,
pm=(_mu_ras*u.mas/u.yr, _mu_decs*u.mas/u.yr),
rv=_vrs*u.km/u.s,
galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
In [ ]:
from gary.dynamics.orbit import combine as combine_orbit
from gary.dynamics.core import combine as combine_ps
In [ ]:
pot = op.load_potential('static_mw')
In [ ]:
w0 = gd.CartesianPhaseSpacePosition(pos=xyz, vel=vxyz)
mean_orbit = pot.integrate_orbit(w0, dt=-0.5, nsteps=12000)
_w0 = gd.CartesianPhaseSpacePosition(pos=samples_xyz, vel=samples_vxyz)
all_w0 = combine_ps((w0,_w0))
orbit = pot.integrate_orbit(all_w0, dt=-0.5, nsteps=12000)
In [ ]:
pers = [orbit[:,i].pericenter().value for i in range(nsamples)] * u.kpc
apos = [orbit[:,i].apocenter().value for i in range(nsamples)] * u.kpc
eccs = (apos - pers) / (apos + pers)
In [ ]:
pers_xyz = np.zeros((3,len(pers)))
pers_xyz[0] = pers.value
pers_xyz = pers_xyz*u.kpc
mx = pot.mass_enclosed(pers_xyz)
pers_xyz = np.zeros((3,len(pers)))
pers_xyz[2] = pers.value
pers_xyz = pers_xyz*u.kpc
mz = pot.mass_enclosed(pers_xyz)
rtide_x = pers.to(u.pc) * (cluster['M'] / (3*mx))**(1/3.)
rtide_z = pers.to(u.pc) * (cluster['M'] / (3*mz))**(1/3.)
peri_rtide = np.mean(np.vstack((rtide_x, rtide_z)).value*rtide_z.unit, axis=0)
core_radius = cluster['Rc']*u.pc
In [ ]:
pl.hist(core_radius / peri_rtide, bins=np.linspace(0,1.,32));
pl.xlabel(r"$r_c / \min(r_t)$")
In [ ]:
pl.hist(eccs, bins=np.linspace(0,1,32))
pl.xlabel("Eccentricity")
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(10,5),sharey=True)
axes[0].hist(pers.value, bins=np.linspace(0,8,32))
axes[0].set_xlabel("Pericenter [kpc]")
axes[1].hist(apos.value, bins=np.linspace(5,20,32))
axes[1].set_xlabel("Apocenter [kpc]")
In [ ]:
period = np.abs(orbit[:,0].estimate_period())
print("orbital period: ", period)
In [ ]:
fig = orbit[:,0].plot()
In [ ]:
c_back,v_back = orbit[:100].to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr, # orbit -50 Myr
galactocentric_frame=galactocentric_frame)
orbit_forw = pot.integrate_orbit(all_w0, dt=0.5, nsteps=100)
c_forw,v_forw = orbit_forw.to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr, # orbit +50 Myr
galactocentric_frame=galactocentric_frame)
In [ ]:
nplot = 256
orbit_alpha = 0.05
back_style = dict(linestyle='-', marker=None, color='k', alpha=orbit_alpha)
forw_style = dict(linestyle='-', marker=None, color='#2166AC', alpha=orbit_alpha)
In [ ]:
ps1_select = mpl.patches.Rectangle((220,-24), width=20, height=7, alpha=0.1)
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(c_back.ra.degree[:,:nplot], c_back.dec.degree[:,:nplot], **back_style)
ax.plot(c_forw.ra.degree[:,:nplot], c_forw.dec.degree[:,:nplot], **forw_style)
ax.set_xlim(ps1_ra_bounds[1].degree+1, ps1_ra_bounds[0].degree-1)
ax.set_ylim(ps1_dec_bounds[0].degree-1, ps1_dec_bounds[1].degree+1)
ax.add_patch(ps1_select)
ax.set_aspect('equal')
pl.minorticks_on()
pl.xlabel("RA [deg]")
pl.ylabel("Dec [deg]");
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(c_back.ra.degree[:,:nplot], c_back.distance[:,:nplot], **back_style)
ax.plot(c_forw.ra.degree[:,:nplot], c_forw.distance[:,:nplot], **forw_style)
ax.set_xlabel("RA [deg]")
ax.set_ylabel("Helio. distance [kpc]")
ax.set_xlim(ps1_ra_bounds[1].degree+1, ps1_ra_bounds[0].degree-1)
ax.set_ylim(9,16)
pl.minorticks_on()
pl.grid()
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(c_back.ra.degree[:,:nplot], distance_modulus(c_back.distance[:,:nplot]), **back_style)
ax.plot(c_forw.ra.degree[:,:nplot], distance_modulus(c_forw.distance[:,:nplot]), **forw_style)
ax.set_xlabel("RA [deg]")
ax.set_ylabel("Distance modulus")
ax.set_xlim(ps1_ra_bounds[1].degree+1, ps1_ra_bounds[0].degree-1)
ax.set_ylim(15.3, 15.9)
# ax.xaxis.set_ticks(np.arange(200,330,10))
pl.minorticks_on()
pl.grid()
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(c_back.ra.degree[:,:nplot], v_back[2].to(u.km/u.s).value[:,:nplot], **back_style)
ax.plot(c_forw.ra.degree[:,:nplot], v_forw[2].to(u.km/u.s).value[:,:nplot], **forw_style)
ax.set_xlabel("RA [deg]")
ax.set_ylabel(r"$v_{\rm los}$ [km/s]")
ax.set_xlim(ps1_ra_bounds[1].degree+1, ps1_ra_bounds[0].degree-1)
pl.minorticks_on()
pl.grid()
Just exclude cluster stars
In [ ]:
control_idx = ps1_c.separation(cluster_c) > out_cluster_pad
targets_idx = np.ones_like(control_idx).astype(bool)
In [ ]:
control = ps1[control_idx]
targets = ps1[targets_idx]
len(control), len(targets)
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(ps1['ra'], ps1['dec'], marker=',', ls='none', alpha=0.1, color='k')
ax.plot(c_back[:,0].ra.degree, c_back[:,0].dec.degree, ls='-', marker=None)
ax.plot(c_forw[:,0].ra.degree, c_forw[:,0].dec.degree, ls='--', marker=None)
ax.set_xlim(ps1_ra_bounds[1].value, ps1_ra_bounds[0].value)
ax.set_ylim(ps1_dec_bounds[0].value, ps1_dec_bounds[1].value)
ax.set_aspect('equal')
pl.minorticks_on()
In [ ]:
cluster_idx = ps1_c.separation(cluster_c) < in_cluster_pad
cluster_stars = ps1[cluster_idx]
len(cluster_stars)
In [ ]:
_buffer = 0.5
fig,ax = pl.subplots(1,1,figsize=(8,6))
alpha = 0.05
ax.plot(ps1['ra'], ps1['dec'], marker='.', ls='none', alpha=alpha, color='k')
ax.plot(cluster_stars['ra'], cluster_stars['dec'], marker='.', ls='none', alpha=0.2, color='g')
ax.set_xlim(cluster['ra']+_buffer,cluster['ra']-_buffer)
ax.set_ylim(cluster['dec']-_buffer,cluster['dec']+_buffer)
ax.set_aspect('equal')
pl.minorticks_on()
In [ ]:
# isochrone
isoch = ascii.read(os.path.join(data_path,"ngc5897_iso_ps1.dat"), header_start=13)
# isoch = ascii.read(os.path.join(data_path,"ngc5897_iso_ps1_brani.dat"), header_start=13)
print(isoch.colnames)
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(6,6))
_iso_style = dict(marker='.', ls='none', alpha=0.5)
ax.plot(isoch['gP1'] - isoch['rP1'], isoch['gP1'] + DM, **_iso_style)
ax.set_xlim(-0.5,1.1)
ax.set_ylim(22,13)
ax.set_xlabel("$(g-r)_0$")
ax.set_ylabel("$g_0$")
fig.tight_layout()
In [ ]:
def plot_color_mag(data, col_prefix="dered_", isochrone=True, dm=15.55, isochrone_style=None, **style):
_g = data['{}g'.format(col_prefix)]
_r = data['{}r'.format(col_prefix)]
_i = data['{}i'.format(col_prefix)]
_z = data['{}z'.format(col_prefix)]
# _y = data['{}y'.format(col_prefix)]
if isochrone_style is None:
isochrone_style = dict()
fig,axes = pl.subplots(1,3,figsize=(8,5),sharey=True)
style['linestyle'] = style.get('linestyle', 'none')
style['marker'] = style.get('marker', '.')
style['alpha'] = style.get('alpha', 0.1)
axes[0].plot(_g-_r, _g, **style)
axes[1].plot(_g-_i, _g, **style)
axes[2].plot(_g-_z, _g, **style)
# axes[3].plot(_g-_y, _g, **style)
if isochrone:
axes[0].plot(isoch['gP1'] - isoch['rP1'], isoch['gP1'] + dm,
marker=isochrone_style.get('marker', None), **isochrone_style)
axes[1].plot(isoch['gP1'] - isoch['iP1'], isoch['gP1'] + dm,
marker=isochrone_style.get('marker', None), **isochrone_style)
axes[2].plot(isoch['gP1'] - isoch['zP1'], isoch['gP1'] + dm,
marker=isochrone_style.get('marker', None), **isochrone_style)
# axes[3].plot(isoch['gP1'] - isoch['yP1'], isoch['gP1'] + dm,
# marker=isochrone_style.get('marker', None), **isochrone_style)
axes[0].set_ylim(22,13)
axes[0].set_xlim(-0.5,1.1)
axes[1].set_xlim(-0.75,2)
axes[2].set_xlim(-1,2.)
# axes[3].set_xlim(-1,2.)
axes[0].set_ylabel("$g_0$")
axes[0].set_xlabel("$(g-r)_0$")
axes[1].set_xlabel("$(g-i)_0$")
axes[2].set_xlabel("$(g-z)_0$")
# axes[3].set_xlabel("$(g-y)_0$")
fig.tight_layout()
return fig, axes
In [ ]:
fig,axes = plot_color_mag(cluster_stars, dm=DM, alpha=0.2)
fig.suptitle("NGC 5897", fontsize=20, y=1.02)
for i,name in enumerate(['g-r','g-i','g-z']):
xy = (min(color_lim[name]),min(g_lim))
w = np.ptp(color_lim[name])
h = np.ptp(g_lim)
r = pl.Rectangle(xy, width=w, height=h, facecolor='g', alpha=0.3)
axes[i].add_patch(r)
In [ ]:
# np.save(os.path.join(data_path, "ngc5897.npy"), cluster_stars)
In [ ]:
def derp(d, colname_pattern="dered_{band}", buffer=0., dm=0):
cp = colname_pattern
b = buffer
g_lim = (13.5, 17.1)
clim = color_lim.copy()
clim['g-r'] = (0,0.8)
g = d[cp.format(band='g')]
gr = g-d[cp.format(band='r')]
gi = g-d[cp.format(band='i')]
gz = g-d[cp.format(band='z')]
ix = (((g+dm) > (g_lim[0]-b)) & ((g+dm) < (g_lim[1]+b)) &
(gr > (clim['g-r'][0]-b)) & (gr < (clim['g-r'][1]+b)))
return d[ix]
glue_ps1 = derp(ps1)
from astropy.table import Table
ascii.write(Table(glue_ps1), "/Users/adrian/projects/globber/data/ngc5897/PS1_bright_glue.txt")
# np.save("/Users/adrian/projects/globber/data/ngc5897/PS1_bright_glue.npy", glue_ps1)
In [ ]:
fig,axes = plot_color_mag(control, dm=DM, alpha=0.05, marker=',', isochrone=False)
fig.suptitle("Control field", fontsize=20, y=1.02)
In [ ]:
fig,axes = plot_color_mag(targets, dm=DM, alpha=0.2, marker=',', isochrone=False)
fig.suptitle("Cluster field", fontsize=20, y=1.02)
In [ ]:
def cut_func(d, colname_pattern="dered_{band}", buffer=0., dm=0):
cp = colname_pattern
b = buffer
g = d[cp.format(band='g')]
gr = g-d[cp.format(band='r')]
gi = g-d[cp.format(band='i')]
gz = g-d[cp.format(band='z')]
ix = (((g+dm) > (g_lim[0]-b)) & ((g+dm) < (g_lim[1]+b)) &
(gr > (color_lim['g-r'][0]-b)) & (gr < (color_lim['g-r'][1]+b)) &
(gi > (color_lim['g-i'][0]-b)) & (gi < (color_lim['g-i'][1]+b)) &
(gz > (color_lim['g-z'][0]-b)) & (gz < (color_lim['g-z'][1]+b)))
return d[ix]
In [ ]:
def draw_lines(axes):
# Just MSTO
# for ax in axes.flat:
# ax.axhline(17.5, color='r')
# ax.axhline(20.5, color='r')
# axes[0].axvline(-0.1, color='r')
# axes[0].axvline(0.6, color='r')
# axes[1].axvline(-0.2, color='r')
# axes[1].axvline(0.9, color='r')
# axes[2].axvline(-0.2, color='r')
# axes[2].axvline(0.9, color='r')
# diagonal lines
# gr_line = lambda x: 10.*x + 15.3
# gi_line = lambda x: 7.*x + 15.25
# gz_line = lambda x: 4.*x + 16.75
# _gr = np.linspace(-0.4,0.8)
# _gi = np.linspace(-0.6,1.5)
# _gz = np.linspace(-0.8,1.5)
# axes[0].plot(_gr, gr_line(_gr), marker=None, color='r')
# axes[1].plot(_gi, gi_line(_gi), marker=None, color='r')
# axes[2].plot(_gz, gz_line(_gz), marker=None, color='r')
# axes[0].plot(_gr, gr_line(_gr)+5, marker=None, color='r')
# axes[1].plot(_gi, gi_line(_gi)+5, marker=None, color='r')
# axes[2].plot(_gz, gz_line(_gz)+4, marker=None, color='r')
# WITH BHB
# for ax in axes.flat:
# ax.axhline(15.5, color='r')
# ax.axhline(20.5, color='r')
# axes[0].axvline(-0.4, color='r')
# axes[0].axvline(0.6, color='r')
# axes[1].axvline(-0.65, color='r')
# axes[1].axvline(0.9, color='r')
# axes[2].axvline(-0.9, color='r')
# axes[2].axvline(0.9, color='r')
# remove the disk instead of selecting the cluster
for ax in axes.flat:
ax.axhline(14., color='r', linestyle='--')
ax.axhline(21.5, color='r', linestyle='--')
ax.axhline(18., color='r')
axes[0].axvline(0.3, color='r')
axes[0].axvline(0.42, color='r')
axes[1].axvline(0.38, color='r')
axes[1].axvline(0.55, color='r')
axes[2].axvline(0.4, color='r')
axes[2].axvline(0.6, color='r')
fig,axes = plot_color_mag(cluster_stars, dm=DM, alpha=0.2)
fig.suptitle("NGC 5897", fontsize=20, y=1.02)
draw_lines(axes)
fig,axes = plot_color_mag(targets, dm=DM, alpha=0.02, marker=',', isochrone=False)
fig.suptitle("Cluster field", fontsize=20, y=1.02)
draw_lines(axes)
In [ ]:
len(cluster_stars)
In [ ]:
cut_targets = cut_func(targets)
cut_isochrone = cut_func(isoch, colname_pattern='{band}P1', dm=DM, buffer=0.05)
cut_control = cut_func(control, buffer=0.05)
cut_cluster = cut_func(cluster_stars, buffer=0.05)
len(cut_targets), len(cut_control)
Hacks below
In [ ]:
def sky_cut(data):
# idx = ((data['ra'] < cluster_c.ra.degree) & (data['ra'] > 225))
idx = ((data['ra'] < cluster_c.ra.degree) & (data['ra'] > 228.4) &
(data['dec'] > cluster_c.dec.degree) & (data['dec'] < -20.5))
return data[idx]
def data_to_X_cov(data):
X = np.vstack([data['dered_{}'.format(band)] for band in 'griz']).T
Xerr = np.vstack([data['{}Err'.format(band)] for band in 'griz']).T
# mixing matrix W
W = np.array([[1, 0, 0, 0], # g magnitude
[1, -1, 0, 0], # g-r color
[1, 0, -1, 0], # g-i color
[1, 0, 0, -1]]) # g-z color
X = np.dot(X, W.T)
# compute error covariance from mixing matrix
Xcov = np.zeros(Xerr.shape + Xerr.shape[-1:])
Xcov[:, range(Xerr.shape[1]), range(Xerr.shape[1])] = Xerr ** 2
# each covariance C = WCW^T
# best way to do this is with a tensor dot-product
Xcov = np.tensordot(np.dot(Xcov, W.T), W, (-2, -1))
return X, Xcov
In [ ]:
clusterX,clusterCov = data_to_X_cov(sky_cut(cluster_stars)[::10])
clusterX.shape
In [ ]:
controlX,controlCov = data_to_X_cov(sky_cut(control))
controlX.shape
In [ ]:
targetsX,targetsCov = data_to_X_cov(sky_cut(targets)[::10])
targetsX.shape
In [ ]:
def kde_with_cov(X, Xcov, Y, Ycov, smooth=0.):
if Xcov.ndim != 4:
Xcov = Xcov[:,np.newaxis,:,:]
# if smooth > 0.:
# H = np.diag(np.ones(Xcov.shape[-1])*smooth)**2
# V = Xcov + H
# else:
V = Xcov + Ycov
LL = log_multivariate_gaussian(X[:,np.newaxis,:], Y, V)
return logsumexp(LL, axis=-1) # NOTE: could also max here
In [ ]:
%timeit LL = kde_with_cov(targetsX, targetsCov, clusterX, clusterCov)
In [ ]:
LL.shape
In [ ]:
this_targets = sky_cut(targets)[::10]
step = 0.05
ra_bins = np.arange(this_targets['ra'].min(), this_targets['ra'].max(), step=step)
dec_bins = np.arange(this_targets['dec'].min(), this_targets['dec'].max(), step=step)
weights = np.exp(LL)
good_idx = np.isfinite(weights)
good_targets = this_targets[good_idx]
H,ra_edges,dec_edges = np.histogram2d(good_targets['ra'], good_targets['dec'],
bins=(ra_bins,dec_bins), weights=weights[good_idx],
normed=True)
# H[H > 1E5] = 1E5
# H_sigma = 1
# H = gaussian_filter(H, sigma=H_sigma)
In [ ]:
pl.imshow(H.T, interpolation='nearest', origin='bottom', cmap='Greys_r', vmin=0.)
In [ ]:
In [ ]:
idx = np.where(LL == LL.max())
In [ ]:
idx
In [ ]:
targetsX[idx[0]], clusterX[idx[1]]
In [ ]:
from sklearn.neighbors import KernelDensity
In [ ]:
features = []
features.append(cut_isochrone['gP1'] + DM)
features.append(cut_isochrone['gP1'] - cut_isochrone['rP1'])
features.append(cut_isochrone['gP1'] - cut_isochrone['iP1'])
features.append(cut_isochrone['gP1'] - cut_isochrone['zP1'])
isoX = np.vstack(features).T
def X_to_grid(X,size=64,custom_limits=None):
nfeatures = X.shape[-1]
grids = []
for i in range(nfeatures):
if custom_limits is None:
minX,maxX = (X[:,i].min(), X[:,i].max())
span = maxX-minX
lims = [minX-0.05*span, maxX+0.05*span]
else:
lims = custom_limits[i]
grid = np.linspace(lims[0], lims[1], size)
grids.append(grid)
return np.meshgrid(*grids)
In [ ]:
# isochrone_kde = KernelDensity(bandwidth=0.05, kernel='gaussian')
isochrone_kde = KernelDensity(bandwidth=0.1, kernel='epanechnikov')
isochrone_kde.fit(isoX)
In [ ]:
Xgrid = X_to_grid(isoX, size=32, custom_limits=[g_lim,color_lim['g-r'],color_lim['g-i'],color_lim['g-z']])
grid_shape = Xgrid[0].shape
Xgrid = np.vstack([X.ravel() for X in Xgrid]).T
log_dens = isochrone_kde.score_samples(Xgrid)
log_dens = np.exp(log_dens).reshape(grid_shape).sum(axis=3).sum(axis=2) / (grid_shape[3]*grid_shape[2])
In [ ]:
pl.figure(figsize=(8,6))
pl.pcolormesh(Xgrid[:,0].reshape(grid_shape)[...,0,0],
Xgrid[:,1].reshape(grid_shape)[...,0,0],
log_dens, cmap='viridis')
pl.plot(cut_isochrone['gP1']+DM, cut_isochrone['gP1']-cut_isochrone['rP1'],
marker='o', alpha=1., linestyle='none', color='r')
# pl.plot(cluster_stars['dered_g']-cluster_stars['dered_r'], cluster_stars['dered_g'],
# marker='.', alpha=0.3, linestyle='none', color='w')
pl.xlim(g_lim[::-1])
pl.ylim(color_lim['g-r'])
pl.colorbar()
pl.gca().set_aspect('equal')
In [ ]:
def kde_with_cov(Xdata, Xcov, Xiso, smooth=None):
if smooth is None:
smooth = isochrone_kde.bandwidth
if Xcov.ndim != 4:
Xcov = Xcov[:,np.newaxis,:,:]
H = np.diag(np.ones(Xcov.shape[-1])*smooth)**2
V = Xcov + H
LL = log_multivariate_gaussian(Xdata[:,np.newaxis,:], Xiso, V)
return logsumexp(LL, axis=-1)
In [ ]:
from astroML.density_estimation import XDGMM
In [ ]:
def data_to_X_cov(data):
X = np.vstack([data['dered_{}'.format(band)] for band in 'griz']).T
Xerr = np.vstack([data['{}Err'.format(band)] for band in 'griz']).T
# mixing matrix W
W = np.array([[1, 0, 0, 0], # g magnitude
[1, -1, 0, 0], # g-r color
[1, 0, -1, 0], # g-i color
[1, 0, 0, -1]]) # g-z color
X = np.dot(X, W.T)
# compute error covariance from mixing matrix
Xcov = np.zeros(Xerr.shape + Xerr.shape[-1:])
Xcov[:, range(Xerr.shape[1]), range(Xerr.shape[1])] = Xerr ** 2
# each covariance C = WCW^T
# best way to do this is with a tensor dot-product
Xcov = np.tensordot(np.dot(Xcov, W.T), W, (-2, -1))
return X, Xcov
In [ ]:
control_X,control_cov = data_to_X_cov(cut_control)
In [ ]:
control_X.shape
TODO: on the cluster, fit an XDGMM with 32 clusters
In [ ]:
# pickle this thing! xd_clf
with open("/Users/adrian/projects/globber/data/xd_control_clf.pickle", "rb") as f:
xd_clf = pickle.load(f)
In [ ]:
# n_clusters = 16
# n_iter = 100
# xd_clf = XDGMM(n_clusters, n_iter=n_iter, tol=1E-4, verbose=True)
# xd_clf.fit(control_X[::20], control_cov[::20])
# # pickle this thing! xd_clf
# with open("/Users/adrian/projects/globber/data/xd_control_clf.pickle", "wb") as f:
# pickle.dump(xd_clf, f)
In [ ]:
sampleX = xd_clf.sample(size=100000)
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(10,5),sharex=True,sharey=True)
axes[0].plot(control_X[::10,1], control_X[::10,0], ls='none', marker=',', alpha=1.)
axes[0].set_title("Data")
axes[1].plot(sampleX[:,1], sampleX[:,0], ls='none', marker=',')
axes[1].set_title("XD")
axes[0].set_xlim(color_lim['g-r'])
axes[0].set_ylim(g_lim)
axes[0].set_ylim(axes[0].get_ylim()[::-1])
Now evaluate with target field
In [ ]:
target_X,target_cov = data_to_X_cov(cut_targets)
In [ ]:
# log_prob_isochrone = isochrone_kde.score_samples(target_X)
log_prob_isochrone = kde_with_cov(target_X, target_cov, isoX, smooth=isochrone_kde.bandwidth)
In [ ]:
log_prob_control = xd_clf.logprob_a(target_X, target_cov)
In [ ]:
log_prob_control.shape
In [ ]:
log_weights = log_prob_isochrone - logsumexp(log_prob_control, axis=-1)
weights = np.exp(log_weights)
In [ ]:
np.isinf(log_weights).sum()
In [ ]:
log_weights.max()
In [ ]:
np.isfinite(log_weights).sum()
In [ ]:
pl.hist(log_weights[np.isfinite(log_weights)], bins=np.linspace(-5,5,128));
In [ ]:
vmin, vmax = (-2, 0)
In [ ]:
pl.figure(figsize=(8,6))
_gr = cut_targets['dered_g']-cut_targets['dered_r']
_g = cut_targets['dered_g']
pl.scatter(_gr[::10], _g[::10], c=log_weights[::10],
marker='.', alpha=0.5, cmap='viridis', vmin=vmin, vmax=vmax)
pl.xlim(-0.5,1.1)
pl.ylim(g_lim[::-1])
pl.colorbar()
In [ ]:
pl.figure(figsize=(8,6))
_ra = cut_targets['ra']
_dec = cut_targets['dec']
pl.scatter(_ra, _dec, c=log_weights, s=4,
marker='o', alpha=0.25, cmap='viridis', vmin=vmin, vmax=vmax)
pl.xlim(233,226)
pl.ylim(-23,-18.5)
pl.gca().set_aspect('equal')
pl.colorbar()
Make filtered image / weighted 2D histogram
In [ ]:
step = 0.08
ra_bins = np.arange(ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree, step=step)
dec_bins = np.arange(ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree, step=step)
In [ ]:
good_idx = np.isfinite(weights)
good_targets = cut_targets[good_idx]
H,ra_edges,dec_edges = np.histogram2d(good_targets['ra'], good_targets['dec'],
bins=(ra_bins,dec_bins), weights=weights[good_idx],
normed=True)
H[H > 1E5] = 1E5
H_sigma = 1
H = gaussian_filter(H, sigma=H_sigma)
In [ ]:
H.max()
In [ ]:
logH = np.log10(H)
fig,axes = pl.subplots(1,3,figsize=(15,5))
axes[0].hist(logH[np.isfinite(logH)], bins=np.linspace(-5,5,64));
axes[0].set_yscale('log')
axes[1].hist(H[np.isfinite(H)], bins=np.linspace(0,0.1,64))
axes[1].set_yscale('log')
In [ ]:
extent = [ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree,
ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree]
fig,axes = pl.subplots(1,3,figsize=(15,5),sharex=True,sharey=True)
axes[0].set_title("log")
axes[0].imshow(logH.T, interpolation='nearest', cmap='Greys',
extent=extent, origin='bottom', vmin=-3, vmax=0)
axes[1].set_title("linear")
axes[1].imshow(H.T, interpolation='nearest', cmap='Greys',
extent=extent, origin='bottom', vmin=0., vmax=0.015)
axes[2].set_title("sqrt")
axes[2].imshow(np.sqrt(H.T), interpolation='nearest', cmap='Greys',
extent=extent, origin='bottom', vmin=0., vmax=np.sqrt(0.015))
for i in range(len(axes)):
axes[i].set_aspect('equal')
# axes[0].set_xlim(235, 225)
axes[0].set_xlim(232, 227)
In [ ]:
def bg_density(data, control, smooth=0.02):
# data
d_g = data['dered_g']
d_gr = data['dered_g'] - data['dered_r']
d_gi = data['dered_g'] - data['dered_i']
# d_gz = data['dered_g'] - data['dered_z']
d_g_var = data['gErr']**2
d_gr_var = d_g_var + data['rErr']**2
d_gi_var = d_g_var + data['iErr']**2
# d_gz_var = d_g_var + data['zErr']**2
# control
c_g = control['dered_g']
c_gr = control['dered_g'] - control['dered_r']
c_gi = control['dered_g'] - control['dered_i']
# c_gz = control['dered_g'] - control['dered_z']
c_g_var = control['gErr']**2
c_gr_var = c_g_var + control['rErr']**2
c_gi_var = c_g_var + control['iErr']**2
# c_gz_var = c_g_var + control['zErr']**2
s_var = smooth**2
const = lambda x,y: 0.5*np.log(2*np.pi*(x + y + s_var))
dist = -0.5 * (d_g[None] - c_g[:,None])**2 / (d_g_var[None] + c_g_var[:,None] + s_var) - const(d_g_var[None], c_g_var[:,None])
dist += -0.5 * (d_gr[None] - c_gr[:,None])**2 / (d_gr_var[None] + c_gr_var[:,None] + s_var) - const(d_gr_var[None], c_gr_var[:,None])
dist += -0.5 * (d_gi[None] - c_gi[:,None])**2 / (d_gi_var[None] + c_gi_var[:,None] + s_var) - const(d_gi_var[None], c_gi_var[:,None])
# dist += -0.5 * (d_gz[None] - c_gz[:,None])**2 / (d_gz_var[None] + c_gz_var[:,None] + s_var) - const(d_gz_var[None], c_gz_var[:,None])
return logsumexp(dist, axis=0) - np.log(len(c_g))
In [ ]:
targets_bg = bg_density(cut_targ, cut_cont[::50], smooth=0.04) # smooth set arbitrarily
In [ ]:
def matched_filter_isochrone(data, isochrone, dm, smooth=0.02):
i_g = isochrone['gP1'] + dm
i_gr = isochrone['gP1'] - isochrone['rP1']
i_gi = isochrone['gP1'] - isochrone['iP1']
# i_gz = isochrone['gP1'] - isochrone['zP1']
d_g = data['dered_g']
d_gr = data['dered_g'] - data['dered_r']
d_gi = data['dered_g'] - data['dered_i']
# d_gz = data['dered_g'] - data['dered_z']
d_g_var = data['gErr']**2
d_gr_var = d_g_var + data['rErr']**2
d_gi_var = d_g_var + data['iErr']**2
# d_gz_var = d_g_var + data['zErr']**2
s_var = smooth**2
const = lambda x: 0.5*np.log(2*np.pi*(x + s_var))
dist = -0.5 * (d_g[None] - i_g[:,None])**2 / (s_var + d_g_var[None]) - const(d_g_var[None])
dist += -0.5 * (d_gr[None] - i_gr[:,None])**2 / (s_var + d_gr_var[None]) - const(d_gr_var[None])
dist += -0.5 * (d_gi[None] - i_gi[:,None])**2 / (s_var + d_gi_var[None]) - const(d_gi_var[None])
# dist += -0.5 * (d_gz[None] - i_gz[:,None])**2 / (s_var + d_gz_var[None]) - const(d_gz_var[None])
return logsumexp(dist, axis=0) - np.log(len(isochrone))
def matched_filter_empirical(data, cluster, dm, smooth=0.02):
# data
d_g = data['dered_g']
d_gr = data['dered_g'] - data['dered_r']
d_gi = data['dered_g'] - data['dered_i']
# d_gz = data['dered_g'] - data['dered_z']
d_g_var = data['gErr']**2
d_gr_var = d_g_var + data['rErr']**2
d_gi_var = d_g_var + data['iErr']**2
# d_gz_var = d_g_var + data['zErr']**2
# control
c_g = cluster['dered_g']
c_gr = cluster['dered_g'] - cluster['dered_r']
c_gi = cluster['dered_g'] - cluster['dered_i']
# c_gz = control['dered_g'] - control['dered_z']
c_g_var = cluster['gErr']**2
c_gr_var = c_g_var + cluster['rErr']**2
c_gi_var = c_g_var + cluster['iErr']**2
# c_gz_var = c_g_var + control['zErr']**2
s_var = smooth**2
const = lambda x,y: 0.5*np.log(2*np.pi*(x + y + s_var))
dist = -0.5 * (d_g[None] - c_g[:,None])**2 / (d_g_var[None] + c_g_var[:,None] + s_var) - const(d_g_var[None], c_g_var[:,None])
dist += -0.5 * (d_gr[None] - c_gr[:,None])**2 / (d_gr_var[None] + c_gr_var[:,None] + s_var) - const(d_gr_var[None], c_gr_var[:,None])
dist += -0.5 * (d_gi[None] - c_gi[:,None])**2 / (d_gi_var[None] + c_gi_var[:,None] + s_var) - const(d_gi_var[None], c_gi_var[:,None])
# dist += -0.5 * (d_gz[None] - c_gz[:,None])**2 / (d_gz_var[None] + c_gz_var[:,None] + s_var) - const(d_gz_var[None], c_gz_var[:,None])
return logsumexp(dist, axis=0) - np.log(len(c_g))
In [ ]:
# filter_smooth = 0.1
# targets_mf = matched_filter_isochrone(cut_targ, cut_isoc,
# dm=dm, smooth=filter_smooth)
targets_mf = matched_filter_empirical(cut_targ, cluster_stars[::2],
dm=dm, smooth=0.05)
In [ ]:
log_weights = targets_mf - targets_bg
log_weights_sort_idx = np.argsort(log_weights)[::-1]
weights = np.exp(log_weights)
In [ ]:
pl.hist(weights, bins=np.logspace(-1,2,64))
pl.xscale('log')
In [ ]:
fig,axes = plot_color_mag(cluster_stars, dm=dm, alpha=0.2)
for j,band in enumerate('riz'):
ix = log_weights_sort_idx[:100]
axes[j].plot(cut_targ[ix]['dered_g']-cut_targ[ix]['dered_{}'.format(band)],
cut_targ[ix]['dered_g'],
marker='.', ls='none')
# for i in range(10):
# for j,band in enumerate('riz'):
# ix = log_weights_sort_idx[i]
# axes[j].plot(cut_targ[ix]['dered_g']-cut_targ[ix]['dered_{}'.format(band)],
# cut_targ[ix]['dered_g'],
# marker='o', color='b')
# ix = log_weights_sort_idx[i+9000]
# axes[j].plot(cut_targ[ix]['dered_g']-cut_targ[ix]['dered_{}'.format(band)],
# cut_targ[ix]['dered_g'],
# marker='o', color='r')
In [ ]:
ra_bins = np.arange(ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree, step=0.1)
dec_bins = np.arange(ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree, step=0.1)
H,ra_edges,dec_edges = np.histogram2d(cut_targ['ra'], cut_targ['dec'], normed=True,
bins=(ra_bins,dec_bins), weights=weights)
_spacing = ra_bins[1]-ra_bins[0]
# H_sigma = (2*u.arcmin).to(u.degree).value / _spacing
H_sigma = 0.6
H = gaussian_filter(H, sigma=H_sigma)
In [ ]:
logH = np.log10(H)
asinhH = np.arcsinh(H)
sqrtH = np.sqrt(H)
In [ ]:
pl.hist(logH[np.isfinite(logH)], bins=np.linspace(-5,5,32));
pl.yscale('log')
In [ ]:
extent = [ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree,
ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree]
pl.imshow(logH.T, interpolation='nearest', cmap='Greys', vmin=-2, vmax=-1.4,
extent=extent, origin='bottom')
pl.xlim(235,225)
In [ ]:
extent = [ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree,
ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree]
pl.imshow(logH.T, interpolation='nearest', cmap='Greys', vmin=-3, vmax=-1.1,
extent=extent, origin='bottom')
pl.xlim(235,225)
In [ ]:
pl.hist(sqrtH.ravel(), bins=np.linspace(0,0.5,32));
In [ ]:
extent = [ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree,
ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree]
pl.imshow(sqrtH.T, interpolation='nearest', cmap='Greys', vmin=0.1, vmax=0.22,
extent=extent, origin='bottom')
pl.xlim(235,225)
In [ ]:
pl.hist(asinhH.ravel(), bins=np.linspace(0,0.1,32));
pl.yscale('log')
In [ ]:
extent = [ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree,
ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree]
pl.imshow(asinhH.T, interpolation='nearest', cmap='Greys', vmin=0.005, vmax=0.04,
extent=extent, origin='bottom')
pl.xlim(235,225)
In [ ]:
pl.hist(H.ravel(), bins=np.linspace(0,0.1,32));
pl.yscale('log')
In [ ]:
extent = [ps1_ra_bounds[0].degree, ps1_ra_bounds[1].degree,
ps1_dec_bounds[0].degree, ps1_dec_bounds[1].degree]
pl.imshow(H.T, interpolation='nearest', cmap='Greys', vmin=0.015, vmax=0.045,
extent=extent, origin='bottom')
pl.xlim(235,225)
In [ ]:
from gary.dynamics.mockstream import fardal_stream, streakline_stream
import gary.integrate as gi
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
for i in range(32):
_plot_filename = os.path.join(plot_path, "{}_fardal.png".format(i))
if os.path.exists(_plot_filename): continue
prog_orbit = pot.integrate_orbit(all_w0[i], dt=-1., nsteps=4000)
stream = fardal_stream(pot, prog_orbit[::-1], prog_mass=2E5,
release_every=2., Integrator=gi.DOPRI853Integrator)
# stream = streakline_stream(pot, prog_orbit[::-1], prog_mass=2E5,
# release_every=2., Integrator=gi.DOPRI853Integrator)
stream_c,stream_v = stream.to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr,
galactocentric_frame=galactocentric_frame)
ax.plot(stream_c.ra.degree, stream_c.dec.degree, marker='.', alpha=0.25, ls='none')
ax.set_xlim(ps1_ra_bounds[1].degree+1, ps1_ra_bounds[0].degree-1)
ax.set_ylim(ps1_dec_bounds[0].degree-1, ps1_dec_bounds[1].degree+1)
ax.set_aspect('equal')
pl.minorticks_on()
pl.xlabel("RA [deg]")
pl.ylabel("Dec [deg]");
fig.savefig(_plot_filename, dpi=300)
ax.cla()
In [ ]: