In [ ]:

    

from collections import OrderedDict
import astropy.units as u
import astropy.coordinates as coord
import numpy as np
from numpy.polynomial.polynomial import polyval
import matplotlib.pyplot as pl
pl.style.use('apw-notebook')
%matplotlib inline
import h5py
from scipy.ndimage import gaussian_filter
from scipy.stats import scoreatpercentile

from astroML.density_estimation import XDGMM
from astroML.utils import log_multivariate_gaussian
from scipy.misc import logsumexp


    

In [ ]:

    

XCov_filename = "/Users/adrian/projects/globber/data/ngc5897/XCov_med.h5"


    

In [ ]:

    

with h5py.File(XCov_filename, "r") as f:
    ra = f['search']['ra'][:]
    dec = f['search']['dec'][:]
    
    allX = f['search']['X'][:]
    allCov = f['search']['Cov'][:]
    
    nonX = f['control']['X'][:]
    nonCov = f['control']['Cov'][:]
    
    cluX = f['cluster']['X'][:]
    isoX = f['isochrone']['X'][:]
    
    non_ll = f['search']['control_log_likelihood'][:]
    iso_ll = f['search']['isochrone_15.55_log_likelihood'][:]
    
    idx = np.isfinite(non_ll)
    
    ra = ra[idx]
    dec = dec[idx]
    allX = allX[idx]
    allCov = allCov[idx]
    non_ll = non_ll[idx]
    iso_ll = iso_ll[idx]


    

In [ ]:

    

bins = np.linspace(-50, 25, 64)
pl.hist(iso_ll, bins=bins, alpha=0.5)
pl.hist(non_ll, bins=bins, alpha=0.5)
# pl.hist(ll, bins=bins, alpha=0.5)
# pl.axvline(threshold)
pl.yscale('log')


    

In [ ]:

    

def worker(allX, allCov, otherX, otherCov, smooth=None, alpha=0.02):
    V = allCov[:,np.newaxis,:,:] + otherCov

    if smooth is not None:
        H = np.zeros(allCov.shape[-2:]) + smooth**2
        V += H[np.newaxis,np.newaxis]

    ll = log_multivariate_gaussian(allX[:,np.newaxis,:], otherX, V)
    N = ll.shape[1]
    bg = N*alpha/(1-alpha) * np.ones((ll.shape[0],1))
    ll = np.hstack((ll,bg))
    ll = logsumexp(ll, axis=-1) - np.log(N*(1-alpha))
    return ll


    

In [ ]:

    

bad_idx, = np.where(non_ll < -10)

bad_ra = ra[bad_idx]
bad_dec = dec[bad_idx]

bad_idx = bad_idx[(bad_ra > ps1['ra'].min()) & (bad_ra < ps1['ra'].max()) &
                  (bad_dec > ps1['dec'].min()) & (bad_dec < ps1['dec'].max())]


    

In [ ]:

    

fig,axes = pl.subplots(1,3,figsize=(12,6),sharex=True,sharey=True)
for i,j in enumerate([1,2,3]):
    axes[i].plot(nonX[::10,j], nonX[::10,0], ls='none', marker=',')

axes[0].set_xlim(0,1.)
axes[0].set_ylim(21,14)


    

In [ ]:

    

fig,axes = pl.subplots(1,3,figsize=(12,6),sharex=True,sharey=True)

# derp_idx = allX[bad_idx,1] < 0.1
for i,j in enumerate([1,2,3]):
    axes[i].plot(allX[bad_idx,j], allX[bad_idx,0], ls='none')

for i,j in enumerate([1,2,3]):
    axes[i].plot(allX[bad_idx[0],j], allX[bad_idx[0],0], ls='none', color='r', marker='o')

axes[0].set_xlim(0,1.)
axes[0].set_ylim(21,14)


    

In [ ]:

    

ll = worker(allX[bad_idx], allCov[bad_idx],
            nonX[::1000], nonCov[::1000], smooth=0.05, alpha=0.1E-2)


    

In [ ]:

    

ll


    

In [ ]:

    

ps1 = np.load("/Users/adrian/projects/globber/data/ngc5897/PS1_stars_pv3_dered_sm.npy")


    

In [ ]:

    

ps1_c = coord.SkyCoord(ra=ps1['ra']*u.degree, dec=ps1['dec']*u.degree)


    

In [ ]:

    

yo = 1
bad_c = coord.SkyCoord(ra=ra[bad_idx[yo]]*u.degree, dec=dec[bad_idx[yo]]*u.degree)
ll[yo].max()


    

In [ ]:

    

match_ix,sep,_ = bad_c.match_to_catalog_sky(ps1_c)
print(sep)
print(ps1[match_ix]['i'] - ps1[match_ix]['dered_i'])


    

In [ ]:

    

((ps1['i'] - ps1['dered_i']) > 0.2).sum()


    

In [ ]:

    

pl.hist(ps1['sg_r'], bins=np.linspace(-6.5,0.2,156))
pl.yscale('log')


    

In [ ]: