In [8]:
import numpy as np
import wisps
import wisps.simulations as wispsim
import pandas as pd
import splat
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord
from scipy.interpolate import interp1d
import numba
from tqdm import tqdm
%matplotlib inline
In [9]:
#dfcand=pd.read_csv(wisps.LIBRARIES+'/candidates.csv')
In [10]:
SIMULATED_DIST=wispsim.simulate_spts()
In [11]:
import wisps.simulations.effective_numbers as ef
In [12]:
data=ef.simulation_outputs()
In [13]:
MASSES=SIMULATED_DIST['mass']
In [14]:
NORM = 0.005/ len(MASSES[np.logical_and(MASSES>0.09, MASSES <=0.1)])
In [15]:
def custom_histogram(things, grid):
n=[]
for g in grid:
n.append(len(things[np.logical_and(g<=things, things< g+1)]))
return np.array(n)
In [16]:
nobs0=custom_histogram((SIMULATED_DIST['spts'][0][:,0]), data['spgrid'])
In [17]:
nobs0
Out[17]:
In [18]:
drop_nan=ef.drop_nan
In [19]:
import splat.empirical as spem
In [20]:
SIMULATED_DIST.keys()
Out[20]:
In [21]:
idx0=np.isnan(SIMULATED_DIST['teffs'])[0]
In [22]:
nan_masses=SIMULATED_DIST['mass'][idx0]
nan_ages=(SIMULATED_DIST['ages'][0])[idx0]
In [23]:
import seaborn as sns
In [24]:
cmap=sns.light_palette((260, 75, 60), input="husl", as_cmap=True)
In [25]:
c=plt.scatter(SIMULATED_DIST['mass'],SIMULATED_DIST['teffs'][4], c=SIMULATED_DIST['ages'][4], marker='+', cmap=cmap, alpha=.5)
plt.xlabel('Mass (Msun)', fontsize=18)
plt.ylabel('Teff ', fontsize=18)
c=plt.colorbar()
c.ax.set_title('Age (Gyr)', fontsize=18)
Out[25]:
In [114]:
from matplotlib.ticker import MultipleLocator
import seaborn as sns
ml1 = MultipleLocator(1.)
ml2 = MultipleLocator(10.)
cmap = sns.cubehelix_palette(reverse=True, as_cmap=True)
In [115]:
cands_df=(wisps.datasets['candidates']).reset_index(drop=True)
In [125]:
#(cands_df[cands_df.spt=='T9.0']).dropna()
In [127]:
extra_cols_cands=((wisps.datasets['stars'])[wisps.datasets['stars'].grism_id.isin(cands_df.grism_id)]).reset_index(drop=True)
In [128]:
cands_df[cands_df.grism_id != extra_cols_cands.grism_id]
Out[128]:
In [129]:
(cands_df[cands_df.spt=='T9.0']).dropna()
Out[129]:
In [130]:
for column in ['survey', 'snr2', 'spt', 'F110', 'F160', 'F140']:
cands_df[column]= extra_cols_cands[column]
In [131]:
cands_df=wisps.Annotator.reformat_table(cands_df)
In [132]:
(cands_df[cands_df.spt=='T9.0']).dropna()
Out[132]:
In [133]:
from matplotlib.colors import Normalize
In [134]:
hs=data['hs']
In [135]:
cnorm=Normalize(hs[0], hs[-1])
In [236]:
#bin in orders of 5 spts
def bin_by_spt_bin(sp_types, number):
ranges=[[17, 20], [20, 25], [25, 30], [30, 35], [35, 40]]
numbers=[]
for r in ranges:
idx= np.logical_and((r[0]<=sp_types), (r[1]>sp_types))
numbers.append(np.nansum(number[idx]))
return numbers
def stay_within_limits(row):
#print (row)
flag=True
if row.survey.lower()=='wisp':
if (row.F110>wisps.MAG_LIMITS['wisps']['F110W'][0] or row.F110<wisps.MAG_LIMITS['wisps']['F110W'][1]):
flag=False
if np.isnan(row.F110):
flag=True
if row.survey.lower()=='hst3d':
if (row.F140>wisps.MAG_LIMITS['hst3d']['F140W'][0] or row.F140<wisps.MAG_LIMITS['hst3d']['F140W'][1]):
flag=False
if np.isnan(row.F140):
flag=True
if splat.typeToNum(row.spt)<17.:
flag=False
return flag
In [231]:
cands_df['stay_flag']=True
In [239]:
flags=cands_df.apply(stay_within_limits, axis=1).values
In [244]:
cdf_to_use=cands_df[flags]
In [245]:
cdf_to_use.shape, cands_df.shape
Out[245]:
In [248]:
nobs=custom_histogram(cdf_to_use.spt.apply(wisps.make_spt_number), data['spgrid'])
In [250]:
nobs.
Out[250]:
In [168]:
cdf_to_use.spt.apply(wisps.make_spt_number)
In [199]:
(cands_df[cands_df.spt.apply(wisps.make_spt_number)==39]).grism_id.dropna().shape
Out[199]:
nobs
In [144]:
nobs.shape, (data['spgrid']).shape
Out[144]:
In [145]:
spgrid2=['M7-L0', 'L0-L5', 'L5-T0', 'T0-T5', 'T5-Y0']
In [154]:
for x, y in zip(data['spgrid'],nobs):
print (x, y)
In [152]:
fig, (ax, ax1)=plt.subplots(ncols=2, figsize=(12, 4))
for x, y in zip(data['spgrid'],nobs):
if y>0:
dy=np.sqrt(y)
dyerr=(y/dy)*np.log10(2.7)*0.5
#ax.errorbar(x,y, yerr=np.sqrt(y), label='observations',fmt='o', color='k')
ax.errorbar(x, np.log10(y), yerr=dyerr, label='observations',fmt='o', color='k')
if y==0:
pass
#ax.plot(x, y, linestyle='none', marker=r'$\downarrow$', label='observations', color='k')
for idx in np.arange(len(hs)):
if hs[idx]==100:
ax.step(data['spgrid'], np.log10(data['n'][idx]*(data['vol'].T)[idx]*NORM), where='mid',
label='h={} pc'.format(hs[idx]),color='#FF851B')
ax1.step(spgrid2, np.log10(bin_by_spt_bin(data['spgrid'], (data['n'][idx])*((data['vol'].T)[idx]*NORM))),
where='mid', label='h={} pc'.format(hs[idx]), color='#FF851B')
elif hs[idx]==1000:
ax.step(data['spgrid'], np.log10(data['n'][idx]*(data['vol'].T)[idx]*NORM),
where='mid', label='h={} pc'.format(hs[idx]), color='#FF851B', linestyle='--')
ax1.step(spgrid2, np.log10(bin_by_spt_bin(data['spgrid'], (data['n'][idx])*((data['vol'].T)[idx]*NORM))),
where='mid', label='h={} pc'.format(hs[idx]),
color='#FF851B', linestyle='--')
else:
ax.step(data['spgrid'], np.log10(data['n'][idx]*(data['vol'].T)[idx]*NORM),
where='mid', label='h={} pc'.format(hs[idx]), color=cmap(cnorm(hs[idx])) )
ax1.step(spgrid2, np.log10(bin_by_spt_bin(data['spgrid'], (data['n'][idx])*((data['vol'].T)[idx]*NORM))),
where='mid', label='h={} pc'.format(hs[idx]), color=cmap(cnorm(hs[idx])))
#ax1.errorbar(spgrid2,bin_by_spt_bin(data['spgrid'],nobs),yerr=np.sqrt(bin_by_spt_bin(data['spgrid'],nobs)),
# label='observations',fmt='o', color='k')
ax1.errorbar(spgrid2,np.log10(bin_by_spt_bin(data['spgrid'],nobs)),
label='observations',fmt='o', color='k')
ax1.legend( fontsize = 15., loc=(1.01, 0.))
ax.minorticks_on()
ax1.minorticks_on()
ax.set_ylim([-.5, 2.1])
ax1.set_ylim([-.5, 2.7])
ax1.tick_params(axis='x', which='minor', bottom=False)
ax.set_ylabel('Log N (SpT)', fontsize=18)
ax1.set_ylabel('Log N (SpT)', fontsize=18)
ax.set_xlabel('SpT', fontsize=18)
ax1.set_xlabel('SpT', fontsize=18)
plt.savefig(wisps.OUTPUT_FIGURES+'/oberved_numbers.pdf', bbox_inches='tight')
In [76]:
wisps.MAG_LIMITS
Out[76]:
In [77]:
co_sl_prob=np.array(data['sl_prob'])
co_dist=np.array(data['dists'])
co_f140=np.array(data['appf140s'])
co_snrj=np.array(data['snrjs'])
In [78]:
co_sl_prob[0].shape, co_dist.shape
Out[78]:
In [79]:
rdf=pd.DataFrame()
rdf['mags0']=np.concatenate(co_f140)
rdf['mags1']=np.concatenate(data['appf160s'])
In [80]:
flags0=rdf.mags0.between(18, 22.5).values
flags1=rdf.mags1.between(16, 21.5).values
In [81]:
flag=np.logical_or(flags0, flags1)
In [82]:
flag.shape, flags0.shape, flags1.shape, rdf.shape
Out[82]:
In [83]:
rdf[flag].shape
Out[83]:
In [84]:
cmap_diverge=sns.diverging_palette(118, 33, n=9, as_cmap=True)
In [85]:
fig, ax=plt.subplots(ncols=3, sharey=True, figsize=(12, 4))
p=ax[0].scatter( np.array(data['appf110s']), np.log10(co_snrj), s=1., marker=',', c=co_sl_prob, cmap='viridis', alpha=0.01)
p=ax[1].scatter( np.array(data['appf140s']), np.log10(co_snrj), s=1., marker=',', c=co_sl_prob, cmap='viridis', alpha=0.01)
p=ax[2].scatter( np.array(data['appf160s']), np.log10(co_snrj), s=1., marker=',', c=co_sl_prob, cmap='viridis', alpha=0.01)
for a in ax:
a.set_ylabel('Log SNR-J', fontsize=18)
a.minorticks_on()
ax[0].set_xlabel('F110W', fontsize=18)
ax[1].set_xlabel('F140W', fontsize=18)
ax[2].set_xlabel('F160W', fontsize=18)
plt.tight_layout()
#ax[0].axhline(21.5, c='k', linestyle='--')
#ax[0].axvline(np.log10(3.), c='k', linestyle='--')
In [86]:
data.keys()
Out[86]:
In [87]:
fig, ax=plt.subplots(ncols=2, nrows=2, figsize=(10, 6))
for idx, h in enumerate(hs):
p=ax[0][0].scatter(np.log10(co_snrj)[idx],drop_nan(SIMULATED_DIST['spts'][0][:,0]), s=1.,
marker=',', c=co_sl_prob[idx], cmap='viridis', alpha=0.01)
p=ax[0][1].scatter(np.log10((data['dists'][idx])),drop_nan(SIMULATED_DIST['spts'][0][:,0]), s=1.,
marker=',', c=co_sl_prob[idx], cmap='viridis', alpha=0.01)
p=ax[1][0].scatter((data['appf140s'][idx]),drop_nan(SIMULATED_DIST['spts'][0][:,0]), s=1.,
marker=',', c=co_sl_prob[idx],cmap='viridis', alpha=0.01)
p=ax[1][1].scatter((data['appf160s'][idx]),drop_nan(SIMULATED_DIST['spts'][0][:,0]), s=1.,
marker=',', c=co_sl_prob[idx],cmap='viridis', alpha=0.01)
ax[1][1].axvline(wisps.MAG_LIMITS['hst3d']['F140W'][0], c='k', linestyle='--')
ax[1][0].axvline(wisps.MAG_LIMITS['hst3d']['F160W'][0], c='k', linestyle='--')
ax[0][0].set_ylabel('SpT', fontsize=18)
ax[0][1].set_ylabel('SpT', fontsize=18)
ax[1][0].set_ylabel('SpT', fontsize=18)
ax[1][1].set_ylabel('SpT', fontsize=18)
plt.tight_layout()
ax[0][0].set_xlabel('Log SNR-J', fontsize=18)
ax[0][1].set_xlabel('Log distance', fontsize=18)
Out[87]:
In [88]:
#the selection function should be the number of objects that I select over the number of objects i "SHOULD select"
#not the number of objects i simulated
In [89]:
from tqdm import tqdm
In [90]:
fig, ax=plt.subplots(ncols=3, nrows=2, figsize=(10, 6))
h= ax[0][0].hist(SIMULATED_DIST['mass'], histtype='step', bins='auto', color='#0074D9')
for idx in tqdm(np.arange(len(hs))):
col=cmap(cnorm(hs[idx]))
linestyle='-'
if hs[idx]==100:
col='#FF851B'
if hs[idx]==1000:
col='#FF851B'
linestyle='--'
h=ax[1][0].hist(drop_nan(np.log10((np.array(data['dists'])[idx]))), bins='auto', histtype='step', label='h={} pc'.format(hs[idx]),
color=col, linestyle= linestyle)
h=ax[1][1].hist(drop_nan(np.log10((np.array(data['snrjs'])[idx]))), bins='auto', histtype='step',label='h={} pc'.format(hs[idx]),
color=col, linestyle= linestyle)
h=ax[1][2].hist(drop_nan((np.array(data['appf140s'])[idx])), bins='auto', histtype='step', label='h={} pc'.format(hs[idx]),
color=cmap(cnorm(hs[idx])))
ax[1][0].set_yscale('log')
ax[1][1].set_yscale('log')
ax[1][2].set_yscale('log')
ax[1][2].axvline(21.5)
h=ax[0][2].hist(SIMULATED_DIST['spts'][0][:,0], histtype='step', bins='auto', color='#0074D9')
h=ax[0][1].hist(SIMULATED_DIST['ages'][0], histtype='step', bins='auto', color='#0074D9')
ax[0][0].set_xlabel(r'Mass (M$_\odot$)', fontsize=18)
ax[0][2].set_xlabel('Spectral Type', fontsize=18)
ax[0][1].set_xlabel('Age (Gyr)', fontsize=18)
ax[1][0].set_xlabel('Log Distance (pc)', fontsize=18)
ax[1][1].set_xlabel('Log SNR-J', fontsize=18)
ax[1][2].set_xlabel('F140W', fontsize=18)
for a in np.concatenate(ax):
a.minorticks_on()
a.set_ylabel('N', fontsize=18)
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/simulations_dists.pdf', bbox_inches='tight')
In [91]:
data.keys()
Out[91]:
In [92]:
import splat.empirical as spe
In [93]:
def ryan_lf(spt):
#ryan's luminosity function
J=spe.typeToMag(spt, '2MASS J')[0]
logphi=-0.30 + 0.11*(J-14) + 0.15*(J -14)**2.+ 0.015*(J-14)**3-0.00020*(J-14)**4
return (10**logphi)*(10**-3)
In [94]:
spgrid=data['spgrid']
In [95]:
(data['vol']).shape
Out[95]:
In [96]:
ryan_n=np.array([ryan_lf(spgrid[idx])*(data['vol'].T)[0][idx] for idx in range(len(spgrid)) ])
In [97]:
ryan_n.shape
Out[97]:
In [98]:
hist=custom_histogram(drop_nan(SIMULATED_DIST['spts'][idx][:,0]), data['spgrid'])
In [99]:
VOLUMES=(data['vol'].T)
In [100]:
(data['n']).shape, (data['vol']).shape, wisps.MAG_LIMITS['hst3d']['F140W']
Out[100]:
In [101]:
(data['spgrid']).shape, (np.log10(data['n'][idx]*(data['vol'].T)[idx]*NORM)).shape
Out[101]:
In [102]:
snew=np.arange(17, 39)
In [103]:
snew
Out[103]:
In [104]:
sptidx=[i for i, spt in enumerate(data['spgrid']) if spt in snew]
In [105]:
#idx
In [106]:
sim_vs=pd.DataFrame(index=[splat.typeToNum(x).replace('.0', ' ') for x in (data['spgrid'])[sptidx]])
for hindex, h in enumerate(hs):
sim_vs['volume {}'.format(h)]=np.round(VOLUMES[hindex])[sptidx]
sim_vs['Number ex {}'.format(h)]=np.round(hist*VOLUMES[hindex]*NORM, 1)[sptidx]
sim_vs['Number ex {}'.format(h)]=(sim_vs['Number ex {}'.format(h)]).astype(int)[sptidx]
sim_vs['Number obs']=nobs[sptidx]
In [107]:
sim_vs.loc["Total"] = sim_vs.sum()
In [108]:
sim_vs
Out[108]:
In [109]:
sim_vs=sim_vs.reindex(sorted(sim_vs.columns), axis=1)
In [110]:
len(sim_vs), sim_vs.shape
Out[110]:
In [111]:
sim_vs.to_latex(wisps.LIBRARIES+'/expectations.tex', index=True)
In [112]:
#(wisps.datasets['candidates']).spt.apply(wisps.make_spt_number).values
In [ ]: