In [ ]:
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.units as u
import h5py
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
pl.style.use('apw-notebook')
%matplotlib inline
from scipy import interpolate
from scipy.misc import logsumexp
from astroML.utils import log_multivariate_gaussian
from globber.ngc5897 import r_t, r_c, cluster_c
from globber.ngc5897 import fiducial_DM as DM
In [ ]:
# XCov_filename = "/Users/adrian/projects/globber/data/ngc5897/XCov.h5"
In [ ]:
# with h5py.File(XCov_filename, "r") as f:
# iso_X = f['isochrone']['X'][:]
# clu_X = f['cluster']['X'][:]
In [ ]:
dart_iso_filename = "/Users/adrian/projects/globber/data/ngc5897/dartmouth_iso_ps1.dat"
dart_iso = ascii.read(dart_iso_filename, header_start=8)
# pars_iso_filename = "/Users/adrian/projects/globber/data/ngc5897/ngc5897_iso_ps1.dat"
# pars_iso = ascii.read(pars_iso_filename, header_start=13)
# pars_iso[114:] = pars_iso[114:][::-1]
pars_iso_filename = "/Users/adrian/projects/globber/data/ngc5897/parsec_iso_ps1.dat"
pars_iso = ascii.read(pars_iso_filename, header_start=13)[:114]
In [ ]:
ps1 = np.load("/Users/adrian/projects/globber/data/ngc5897/PS1_stars_pv3_dered_sm.npy")
ps1_c = coord.SkyCoord(ra=ps1['ra']*u.degree, dec=ps1['dec']*u.degree)
cluster = ps1[ps1_c.separation(cluster_c) < 3*r_c]
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(8,8),sharex=True,sharey=True)
axes[0].plot(cluster['dered_g']-cluster['dered_i'], cluster['dered_i'], ls='none', marker='.', alpha=0.5)
axes[0].plot(dart_iso['gP1']-dart_iso['iP1'], dart_iso['iP1']+DM, marker=None, lw=2.)
axes[0].plot(pars_iso['gP1']-pars_iso['iP1'], pars_iso['iP1']+DM, marker=None, lw=2.)
# axes[1].plot(clu_X[:,2], clu_X[:,0], ls='none', alpha=0.5)
axes[1].plot(dart_iso['gP1']-dart_iso['iP1'], dart_iso['iP1']+DM, marker=None, lw=2.)
axes[1].plot(pars_iso['gP1']-pars_iso['iP1'], pars_iso['iP1']+DM, marker=None, lw=2.)
pl.xlim(0,0.7)
pl.ylim(21,17)
In [ ]:
itrp_i = np.linspace(17., 21.5, 128)
itrp_gi = interp(itrp_i)
In [ ]:
pl.figure(figsize=(4,8))
pl.scatter(itrp_gi, itrp_i, alpha=0.5)
pl.xlim(0,1.)
pl.ylim(21,17.5)
In [ ]:
x = np.arange(0,1,0.02)
y = np.arange(17.5,21.,0.04)
shp = (y.size, x.size)
xygrid = np.vstack(list(map(np.ravel,np.meshgrid(x, y)))).T
In [ ]:
xymodel = np.vstack((itrp_gi,itrp_i)).T
xymodel.shape
In [ ]:
h = 0.02
V = np.diag([h]*2)**2
W = np.array([[1, -1], # g-i
[0, 1]]) # i magnitude
# each covariance C = WCW^T
# V = np.einsum('mj,jk->mk', W, V)
# V = np.einsum('lk,mk->ml', W, V)
# V
In [ ]:
ll = log_multivariate_gaussian(xygrid[:,np.newaxis], xymodel[np.newaxis], V)
ll = logsumexp(ll, axis=1)
ll.shape
In [ ]:
pl.figure(figsize=(4,8))
pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp),
np.exp(ll).reshape(shp), cmap='Blues')
pl.xlim(0,1.)
pl.ylim(21,17.5)
In [ ]:
with h5py.File(XCov_filename, "r") as f:
bg_X = f['control']['X'][:]
bg_Cov = f['control']['Cov'][:]
idx = (bg_X[:,0] >= 17.) & (bg_X[:,0] <= 21.5) & (bg_X[:,2] >= -0.1) & (bg_X[:,2] <= 1.1)
bg_X = bg_X[idx]
bg_Cov = bg_Cov[idx]
bg_X = bg_X[::50,[2,0]]
bg_Cov = bg_Cov[::50,[2,0]][:,:,[2,0]]
In [ ]:
pl.figure(figsize=(4,8))
pl.plot(bg_X[:,0], bg_X[:,1], alpha=0.4, marker='.', ls='none')
pl.xlim(0,1.)
pl.ylim(21,17.5)
In [ ]:
bg_h = 0.08
bg_V = np.diag([bg_h]*2)**2
W = np.array([[1, -1], # g-i
[0, 1]]) # i magnitude
# each covariance C = WCW^T
bg_V = np.einsum('mj,jk->mk', W, bg_V)
bg_V = np.einsum('lk,mk->ml', W, bg_V)
bg_V
In [ ]:
bg_X[np.newaxis].shape, xygrid[:,np.newaxis].shape, _V[np.newaxis].shape
In [ ]:
_V = bg_Cov + bg_V[np.newaxis]
bg_ll = log_multivariate_gaussian(bg_X[np.newaxis], xygrid[:,np.newaxis], _V[np.newaxis])
bg_ll = logsumexp(bg_ll, axis=1)
bg_ll.shape
In [ ]:
pl.figure(figsize=(4,8))
pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp),
np.exp(bg_ll).reshape(shp), cmap='Blues')
pl.xlim(0,1.)
pl.ylim(21,17.5)
In [ ]:
pl.figure(figsize=(4,8))
pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp),
np.exp(ll - bg_ll).reshape(shp), cmap='Blues')
pl.xlim(0,1.)
pl.ylim(21,17.5)
In [ ]:
_iso = iso[::-1]
tck, u = interpolate.splprep([_iso['gP1']-_iso['rP1'],_iso['gP1']-_iso['iP1'],_iso['gP1']-_iso['zP1']],
u=_iso['iP1'], k=3, s=1E-4)
u_fine = np.linspace(u.min(), u.max(), 1024)
gr,gi,gz = interpolate.splev(u_fine, tck)
# tck, u = interpolate.splprep([_iso['gP1']-_iso['iP1']],
# u=_iso['iP1'], k=3, s=1E-4)
# u_fine = np.linspace(u.min(), u.max(), 1024)
# gi = interpolate.splev(u_fine, tck)
In [ ]:
pl.figure(figsize=(4,8))
pl.scatter(itrp_gi, itrp_i, alpha=0.5, color='k')
pl.scatter(gz, u_fine+DM, alpha=0.5, color='r')
pl.scatter(gi, u_fine+DM, alpha=0.5, color='g')
pl.scatter(gr, u_fine+DM, alpha=0.5, color='b')
pl.xlim(0,1.)
pl.ylim(21,17.5)
In [ ]:
bg_h = 0.05
bg_V = np.diag([bg_h]*2)**2
W = np.array([[1, -1], # g-i
[0, 1]]) # i magnitude
# each covariance C = WCW^T
bg_V = np.einsum('mj,jk->mk', W, bg_V)
bg_V = np.einsum('lk,mk->ml', W, bg_V)
bg_V
In [ ]:
_Cov = bg_Cov[:10]
_X = bg_X[:10]
_Cov.shape, _X.shape
In [ ]:
_X2 = xygrid
_Cov2 = np.zeros(_X2.shape + (_X2.shape[-1],))
_Cov2.shape, _X2.shape
In [ ]:
# METHOD 1:
_V = _Cov + bg_V[np.newaxis] # uncertainties plus smooth
_ll1 = log_multivariate_gaussian(_X[np.newaxis], _X2[:,np.newaxis], _V[np.newaxis])
_ll1 = logsumexp(_ll1, axis=1)
_ll1.shape
In [ ]:
# METHOD 2:
_V12 = _Cov[np.newaxis] + _Cov2[:,np.newaxis] + bg_V[np.newaxis,np.newaxis] # uncertainties plus smooth
_ll2 = log_multivariate_gaussian(_X[np.newaxis], _X2[:,np.newaxis], _V)
_ll2 = logsumexp(_ll2, axis=1)
_ll2.shape
# print(_V12.shape, _X.shape, _X2.shape, _ll2.shape)
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(8,8),sharex=True,sharey=True)
axes[0].pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp),
np.exp(_ll1).reshape(shp), cmap='Blues')
axes[1].pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp),
np.exp(_ll2).reshape(shp), cmap='Blues')
pl.xlim(0,1.)
pl.ylim(21,17.5)
In [ ]:
def worker(allX, allCov, otherX, otherCov=None, smooth=None):
if otherCov is not None:
V = allCov[np.newaxis] + otherCov[:,np.newaxis]
if smooth is not None:
H = np.eye(allCov.shape[-1]) * smooth**2
V += H[np.newaxis,np.newaxis]
ll = log_multivariate_gaussian(allX[np.newaxis], otherX[:,np.newaxis], V)
else:
V = allCov
if smooth is not None:
H = np.eye(allCov.shape[-1]) * smooth**2
V += H[np.newaxis]
ll = log_multivariate_gaussian(allX[np.newaxis], otherX[:,np.newaxis], V[np.newaxis])
ll = logsumexp(ll, axis=0) # NOTE: could also max here
return ll
In [ ]:
worker(_X, _Cov, _X2).shape
In [ ]:
with h5py.File(XCov_filename, "r") as f:
_X = f['isochrone']['X'][:]
In [ ]:
fig,axes = pl.subplots(1,4,figsize=(15,5), sharey=True)
DERPS = [_iso['gP1']-_iso['rP1'],_iso['gP1']-_iso['iP1'],_iso['gP1']-_iso['zP1'],_iso['rP1']-_iso['zP1']]
for i in range(1,4+1):
axes[i-1].plot(_X[:,i], _X[:,0]+DM, ls='none', alpha=0.4)
axes[i-1].plot(DERPS[i-1], _iso['iP1']+DM, ls='none', alpha=0.4)
axes[0].set_ylim(22,14)
In [ ]: