<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Now a distance cut
Step1:  Stars I actually observed
Step2:  Data for the observed stars
Step3:  Comparison of stars observed with Catalina
Step4:  Issues
Why are some of the PTF RR Lyrae missing from Catalina? Because they are too faint! (R>18)
Why are Catalina stars missing from PTF? More observations, larger selection window.
Step5:  Possible Blaschko stars
Step6:  For Kathryn's proposal
Step7:  Now read in RR Lyrae data, compute prob for each star
<ASSISTANT_TASK:>
Python Code:
d = triand['dh'].data
d_cut = (d > 15) & (d < 21)
triand_dist = triand[d_cut]
c_triand = _c_triand[d_cut]
print(len(triand_dist))
plt.hist(triand_dist['<Vmag>'].data)
Explanation: Now a distance cut:
End of explanation
ptf_triand = ascii.read("/Users/adrian/projects/streams/data/observing/triand.txt")
ptf_c = coord.SkyCoord(ra=ptf_triand['ra']*u.deg, dec=ptf_triand['dec']*u.deg)
print ptf_triand.colnames, len(ptf_triand)
obs_dist = distance(ptf_triand['Vmag'].data)
((obs_dist > 12*u.kpc) & (obs_dist < 25*u.kpc)).sum()
ptf_triand[0]
Explanation: Stars I actually observed
End of explanation
rrlyr_d = np.genfromtxt("/Users/adrian/projects/triand-rrlyrae/data/RRL_ALL.txt", 
                        skiprows=2, dtype=None, names=['l','b','vhel','vgsr','src','ra','dec','name','dist'])
obs_rrlyr = rrlyr_d[rrlyr_d['src'] == 'PTF']
Explanation: Data for the observed stars
End of explanation
fig,ax = plt.subplots(1,1,figsize=(10,8))
# ax.plot(c.galactic.l.degree, c.galactic.b.degree, linestyle='none',
#         marker='o', markersize=4, alpha=0.75) # ALL RR LYRAE
ax.plot(c_triand.galactic.l.degree, c_triand.galactic.b.degree, linestyle='none',
        marker='o', markersize=5, alpha=0.75)
ax.plot(ptf_c.galactic.l.degree, ptf_c.galactic.b.degree, linestyle='none', 
        marker='o', markerfacecolor='none', markeredgewidth=2, markersize=12, alpha=0.75)
ax.plot(obs_rrlyr['l'], obs_rrlyr['b'], linestyle='none', mec='r',
        marker='o', markerfacecolor='none', markeredgewidth=2, markersize=12, alpha=0.75)
# x = np.linspace(-10,40,100)
# x[x < 0] += 360.
# y = np.linspace(30,45,100)
# x,y = map(np.ravel, np.meshgrid(x,y))
# ccc = coord.SkyCoord(ra=x*u.deg,dec=y*u.deg)
# ax.plot(ccc.galactic.l.degree, ccc.galactic.b.degree, linestyle='none')
ax.set_xlim(97,162)
ax.set_ylim(-37,-13)
ax.set_xlabel("$l$ [deg]")
ax.set_ylabel("$b$ [deg]")
Explanation: Comparison of stars observed with Catalina
End of explanation
fig,ax = plt.subplots(1,1,figsize=(10,8))
ax.plot(c_triand.galactic.l.degree, c_triand.galactic.b.degree, linestyle='none',
        marker='o', markersize=4, alpha=0.75)
ax.plot(ptf_c.galactic.l.degree, ptf_c.galactic.b.degree, linestyle='none', 
        marker='o', markerfacecolor='none', markeredgewidth=2, markersize=8, alpha=0.75)
ax.plot(obs_rrlyr['l'], obs_rrlyr['b'], linestyle='none', mec='r',
        marker='o', markerfacecolor='none', markeredgewidth=2, markersize=8, alpha=0.75)
ax.plot(c_triand.galactic.l.degree[10], c_triand.galactic.b.degree[10], linestyle='none',
        marker='o', markersize=25, alpha=0.75)
ax.set_xlim(97,162)
ax.set_ylim(-37,-13)
c_triand.icrs[10]
Explanation: Issues
Why are some of the PTF RR Lyrae missing from Catalina? Because they are too faint! (R>18)
Why are Catalina stars missing from PTF? More observations, larger selection window.
End of explanation
brani = ascii.read("/Users/adrian/projects/triand-rrlyrae/brani_sample/TriAnd.dat")
blaschko = brani[(brani['objectID'] == "13322281016459551106") | (brani['objectID'] == "13879390364114107826")]
for b in blaschko:
    row = ptf_triand[np.argmin(np.sqrt((ptf_triand['ra'] - b['ra'])**2 + (ptf_triand['dec'] - b['dec'])**2))]
    print(row['name'])
    print(coord.SkyCoord(ra=row['ra']*u.deg, dec=row['dec']*u.deg).galactic)
zip(obs_rrlyr['l'], obs_rrlyr['b'])
d = V_to_dist(triand['<Vmag>'].data).to(u.kpc).value
bins = np.arange(1., 60+5, 3)
plt.figure(figsize=(10,8))
n,bins,patches = plt.hist(triand['dh'].data, bins=bins, alpha=0.5, label='Catalina')
for pa in patches:
    if pa.xy[0] < 15. or pa.xy[0] > 40.:
        pa.set_alpha(0.2)
# other_bins = np.arange(0, 15+2., 2.)
# plt.hist(V_to_dist(triand['<Vmag>'].data), bins=other_bins, alpha=0.2, color='k')
# other_bins = np.arange(40, 60., 2.)
# plt.hist(V_to_dist(triand['<Vmag>'].data), bins=other_bins, alpha=0.2, color='k')
plt.hist(V_to_dist(ptf_triand['Vmag'].data), 
         bins=bins, alpha=0.5, label='PTF/MDM')
plt.xlabel("Distance [kpc]")
plt.ylabel("Number")
# plt.ylim(0,35)
plt.legend(fontsize=20)
plt.axvline(18.)
plt.axvline(28.)
Explanation: Possible Blaschko stars:
* R_13322281016459551106
* R_13879390364114107826
End of explanation
import emcee
import triangle
from scipy.misc import logsumexp
((distance(triand['<Vmag>'].data) > (15.*u.kpc)) & (distance(triand['<Vmag>'].data) < (40.*u.kpc))).sum()
!head -n3 /Users/adrian/projects/triand-rrlyrae/data/triand_giants.txt
d = np.loadtxt("/Users/adrian/projects/triand-rrlyrae/data/triand_giants.txt", skiprows=1)
d2 = np.genfromtxt("/Users/adrian/projects/triand-rrlyrae/data/TriAnd_Mgiant.txt", skiprows=2)
plt.plot(d[:,0], d[:,2], linestyle='none')
plt.plot(d2[:,0], d2[:,3], linestyle='none')
ix = (d[:,2] < 100) & (d[:,2] > -50)
ix = np.ones_like(ix).astype(bool)
plt.plot(d[ix,0], d[ix,2], linestyle='none')
plt.plot(d[ix,0], -1*d[ix,0] + 170, marker=None)
plt.xlabel('l [deg]')
plt.ylabel('v_r [km/s]')
plt.figure()
plt.plot(d[ix,0], d[ix,1], linestyle='none')
plt.xlabel('l [deg]')
plt.ylabel('b [deg]')
def ln_normal(x, mu, sigma):
    return -0.5*np.log(2*np.pi) - np.log(sigma) - 0.5*((x-mu)/sigma)**2
# def ln_prior(p):
#     m,b,V = p
    
#     if m > 0. or m < -50:
#         return -np.inf
    
#     if b < 0 or b > 500:
#         return -np.inf
    
#     if V <= 0.:
#         return -np.inf
    
#     return -np.log(V)
# def ln_likelihood(p, l, vr, sigma_vr):
#     m,b,V = p
#     sigma = np.sqrt(sigma_vr**2 + V**2)
#     return ln_normal(vr, m*l + b, sigma)
# mixture model - f_ol is outlier fraction
def ln_prior(p):
    m,b,V,f_ol = p
    
    if m > 0. or m < -50:
        return -np.inf
    
    if b < 0 or b > 500:
        return -np.inf
    
    if V <= 0.:
        return -np.inf
    
    if f_ol > 1. or f_ol < 0.:
        return -np.inf
    
    return -np.log(V)
def likelihood(p, l, vr, sigma_vr):
    m,b,V,f_ol = p
    sigma = np.sqrt(sigma_vr**2 + V**2)
    term1 = ln_normal(vr, m*l + b, sigma)
    term2 = ln_normal(vr, 0., 120.)
    return np.array([term1, term2])
def ln_likelihood(p, *args):
    m,b,V,f_ol = p
    x = likelihood(p, *args)
    
    # coefficients
    b = np.zeros_like(x)
    b[0] = 1-f_ol
    b[1] = f_ol
    
    return logsumexp(x,b=b, axis=0)
    
def ln_posterior(p, *args):
    lnp = ln_prior(p)
    if np.isinf(lnp):
        return -np.inf
    
    return lnp + ln_likelihood(p, *args).sum()
def outlier_prob(p, *args):
    m,b,V,f_ol = p
    p1,p2 = likelihood(p, *args)
    return f_ol*np.exp(p2) / ((1-f_ol)*np.exp(p1) + f_ol*np.exp(p2))
vr_err = 2 # km/s
nwalkers = 32
sampler = emcee.EnsembleSampler(nwalkers=nwalkers, dim=4, lnpostfn=ln_posterior, 
                                args=(d[ix,0],d[ix,2],vr_err))
p0 = np.zeros((nwalkers,sampler.dim))
p0[:,0] = np.random.normal(-1, 0.1, size=nwalkers)
p0[:,1] = np.random.normal(150, 0.1, size=nwalkers)
p0[:,2] = np.random.normal(25, 0.5, size=nwalkers)
p0[:,3] = np.random.normal(0.1, 0.01, size=nwalkers)
for pp in p0:
    lnp = ln_posterior(pp, *sampler.args)
    if not np.isfinite(lnp):
        print("you suck")
pos,prob,state = sampler.run_mcmc(p0, N=100)
sampler.reset()
pos,prob,state = sampler.run_mcmc(pos, N=1000)
fig = triangle.corner(sampler.flatchain,
                      labels=[r'$\mathrm{d}v/\mathrm{d}l$', r'$v_0$', r'$\sigma_v$', r'$f_{\rm halo}$'])
figsize = (12,8)
MAP = sampler.flatchain[sampler.flatlnprobability.argmax()]
pout = outlier_prob(MAP, d[ix,0], d[ix,2], vr_err)
plt.figure(figsize=figsize)
cl = plt.scatter(d[ix,0], d[ix,2], c=(1-pout), s=30, cmap='RdYlGn', vmin=0, vmax=1)
cbar = plt.colorbar(cl)
cbar.set_clim(0,1)
# plt.plot(d[ix,0], d[ix,2], linestyle='none', marker='o', ms=4)
plt.xlabel(r'$l\,[{\rm deg}]$')
plt.ylabel(r'$v_r\,[{\rm km\,s}^{-1}]$')
ls = np.linspace(d[ix,0].min(), d[ix,0].max(), 100)
for i in np.random.randint(len(sampler.flatchain), size=100):
    m,b,V,f_ol = sampler.flatchain[i]
    plt.plot(ls, m*ls+b, color='#555555', alpha=0.1, marker=None)
best_m,best_b,best_V,best_f_ol = MAP
plt.plot(ls, best_m*ls + best_b, color='k', alpha=1, marker=None)
plt.plot(ls, best_m*ls + best_b + best_V, color='k', alpha=1, marker=None, linestyle='--')
plt.plot(ls, best_m*ls + best_b - best_V, color='k', alpha=1, marker=None, linestyle='--')
plt.xlim(ls.max()+2, ls.min()-2)
plt.title("{:.1f}% halo stars".format(best_f_ol*100.))
print(((1-pout) > 0.75).tolist())
print best_m, best_b, best_V
print "MAP velocity dispersion: {:.2f} km/s".format(best_V)
high_p = (1-pout) > 0.8
plt.figure(figsize=figsize)
cl = plt.scatter(d[high_p,0], d[high_p,1], c=d[high_p,2]-d[high_p,2].mean(), s=30, cmap='coolwarm', vmin=-40, vmax=40)
cbar = plt.colorbar(cl)
ax = plt.gca()
ax.set_axis_bgcolor('#555555')
plt.xlim(ls.max()+2,ls.min()-2)
plt.ylim(-50,-10)
plt.xlabel(r'$l\,[{\rm deg}]$')
plt.ylabel(r'$b\,[{\rm deg}]$')
plt.title(r'$P_{\rm TriAnd} > 0.8$', y=1.02)
Explanation: For Kathryn's proposal
End of explanation
rrlyr_d = np.genfromtxt("/Users/adrian/projects/triand-rrlyrae/data/RRL_ALL.txt", skiprows=2, dtype=None)
!cat "/Users/adrian/projects/triand-rrlyrae/data/RRL_ALL.txt"
rrlyr_d = np.genfromtxt("/Users/adrian/projects/triand-rrlyrae/data/RRL_ALL.txt", skiprows=2)
rrlyr_vr_err = 10.
MAP = sampler.flatchain[sampler.flatlnprobability.argmax()]
pout = outlier_prob(MAP, rrlyr_d[:,0], rrlyr_d[:,3], rrlyr_vr_err)
plt.figure(figsize=figsize)
cl = plt.scatter(rrlyr_d[:,0], rrlyr_d[:,1], c=(1-pout), s=30, cmap='RdYlGn', vmin=0, vmax=1)
cbar = plt.colorbar(cl)
cbar.set_clim(0,1)
# plt.plot(d[ix,0], d[ix,2], linestyle='none', marker='o', ms=4)
plt.xlabel(r'$l\,[{\rm deg}]$')
plt.ylabel(r'$b\,[{\rm deg}]$')
plt.xlim(ls.max()+2,ls.min()-2)
plt.ylim(-50,-10)
plt.title("RR Lyrae")
MAP = sampler.flatchain[sampler.flatlnprobability.argmax()]
pout = outlier_prob(MAP, rrlyr_d[:,0], rrlyr_d[:,3], rrlyr_vr_err)
plt.figure(figsize=figsize)
cl = plt.scatter(rrlyr_d[:,0], rrlyr_d[:,3], c=(1-pout), s=30, cmap='RdYlGn', vmin=0, vmax=1)
cbar = plt.colorbar(cl)
cbar.set_clim(0,1)
# plt.plot(d[ix,0], d[ix,2], linestyle='none', marker='o', ms=4)
plt.xlabel(r'$l\,[{\rm deg}]$')
plt.ylabel(r'$v_r\,[{\rm km\,s}^{-1}]$')
ls = np.linspace(d[ix,0].min(), d[ix,0].max(), 100)
best_m,best_b,best_V,best_f_ol = MAP
plt.plot(ls, best_m*ls + best_b, color='k', alpha=1, marker=None)
plt.plot(ls, best_m*ls + best_b + best_V, color='k', alpha=1, marker=None, linestyle='--')
plt.plot(ls, best_m*ls + best_b - best_V, color='k', alpha=1, marker=None, linestyle='--')
plt.xlim(ls.max()+2, ls.min()-2)
plt.title("RR Lyrae")
Explanation: Now read in RR Lyrae data, compute prob for each star
End of explanation
<END_TASK>