In [109]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from collections import OrderedDict as odict
Load the SN light curve fits that were generated in this notebook.
This is an array of "SN fit" objects. An example is printed at the end of this cell:
In [110]:
import sys
import gzip, pickle
if sys.version.startswith('2'):
snFits = pickle.load(gzip.GzipFile('snFits.p.gz'))
else:
snFits = pickle.load(gzip.GzipFile('snFits.p.gz'),
encoding='latin1')
print(len(snFits))
snf = [s for s in snFits.values() if s is not None]
print(len(snf))
snf[0]
Out[110]:
Use list comprehension to extract the redshift (z) and fitted distance modulus (mu) from each of the supernova light curve fits. Each of those becomes its own array of numbers.
Also extract mu_var. Need to do this in a loop since for some light curve fits mu_var doesn't exist (for a reason that I am not sure of right now).
In [111]:
z = np.array([s.z for s in snf])
mu = np.array([s.mu for s in snf])
mu_var = mu.copy() #np.array([s.mu_var for s in snf])
for i,s in enumerate(snf):
try:
mu_var[i] = s.mu_var
except:
mu_var[i] = 999
# Here, we didn't fit z so zz and z will be identical.
zz = z.copy() #np.array([s['inputParams'].get('z') for s in snf])
for i,s in enumerate(snf):
try:
zz[i] = s['inputParams'].get('z')
except:
zz[i] = 999
Filter the data to remove points that are probably way wrong (e.g. too high or too low mu).
TODO: we need to look at some of these bad points and understand why their fits gave incorrect results.
In [112]:
z = z[(mu > 0) & (mu < 19)]
zz = zz[(mu > 0) & (mu < 19)]
mu = mu[(mu > 0) & (mu < 19)]
mu_var = mu_var[(mu > 0) & (mu < 19)]
Plot the Hubble diagram.
In [113]:
#plt.scatter(z, mu)
plt.errorbar(z, mu, yerr=np.sqrt(mu_var), fmt='o', ecolor='r', elinewidth=1, ms=2)
plt.xlabel('z')
plt.xscale('log')
plt.ylabel('mu + const')
plt.xlim((0.2,1.4))
plt.ylim((8,18));
The mus computed above are offset by some constant that we do not know yet. Let's estimate that offset by computing a Hubble diagram for a model cosmology.
In [114]:
from astropy.cosmology import Planck15 as cosmo
mu = np.array([s.mu for s in snf])
x0 = np.array([s.x0 for s in snf])
x0 = x0[(mu > 0) & (mu < 19)]
x1 = np.array([s.x1 for s in snf])
x1 = x1[(mu > 0) & (mu < 19)]
c = np.array([s.c for s in snf])
c = c[(mu > 0) & (mu < 19)]
mu = mu[(mu > 0) & (mu < 19)]
muz = np.array(cosmo.distmod(z))
alpha = 0.14
beta = -3.11
M = -2.5 * np.log10(x0) + alpha * x1 + beta * c - muz
# plt.scatter(z, muz)
# plt.xlabel('z')
# plt.xscale('log')
# plt.xlim((0.2,1.4))
# plt.ylabel('mu')
# plt.ylim((39,45));
Compute the constant offset as the median of the differences between the cosmological mu and the fitted SN mu.
In [115]:
const = np.median(muz-mu)
print(const)
plt.scatter(z,muz-mu)
plt.plot([0,1.4],[const,const],'r')
Out[115]:
In [116]:
plt.errorbar(z, mu+const, yerr=np.sqrt(mu_var), fmt='o', ecolor='r', elinewidth=1, ms=2)
plt.xlabel('z')
plt.xscale('log')
plt.ylabel('mu + const')
plt.xlim((0.2,1.4))
plt.xticks(np.arange(0.2,1.4,0.2))
plt.ylim((39,48));
In [117]:
plt.plot(z, np.sqrt(mu_var), 'o')
Out[117]:
In [118]:
import pandas as pd
In [119]:
df= pd.DataFrame(snFits)
In [120]:
df.head()
Out[120]:
In [121]:
df.transpose()
df=df.transpose()
In [122]:
results= df[['z','mu','mu_var']]
results.head()
Out[122]:
In [123]:
lowz= results.query('z<.4')
In [124]:
_ = plt.hist(lowz['mu_var'], bins=np.arange(0., 15, 0.5), histtype='step', normed=1)
In [125]:
highz= results.query('z>.7')
In [126]:
_ = plt.hist(highz['mu_var'], bins=np.arange(0., 15, 0.5), histtype='step', normed=1)
In [127]:
results['bindex']= results['z']//0.1
In [128]:
results.head()
Out[128]:
In [129]:
results[['z', 'mu_var', 'mu']] = results[['z', 'mu_var', 'mu']].astype(np.float)
In [130]:
dic = dict(mu_var=np.mean, )
In [131]:
mean=results.groupby('bindex').agg(dic)
mean
Out[131]:
In [132]:
#sd= dict(mu_var=np.std, )
In [133]:
results.groupby('bindex').agg(sd)
In [134]:
results['mu_err']=np.sqrt(results['mu_var'])
In [135]:
results.head()
Out[135]:
In [136]:
def clipped_mean(mu_err, outlier_threshold):
return np.mean(mu_err['mu_err<outlier_threshold'])
In [137]:
grouped=results.groupby('bindex')
In [138]:
mu_err_table=grouped.agg(dict(mu_err=['count', np.mean, np.std]))
mu_err_table
Out[138]:
In [139]:
mu_err_table['frequencies'] = mu_err_table.mu_err['count'] / mu_err_table.mu_err['count'].sum()
In [140]:
mu_err_table= mu_err_table.ix[1:]
mu_err_table
Out[140]:
In [162]:
NumSN= 50000
numObjectsPerBin = np.round(mu_err_table.frequencies * NumSN).astype(np.int)
print(numObjectsPerBin)
In [142]:
numObjectsPerBin.ix[1]
Out[142]:
In [143]:
mu_err_table['mu_err'].ix[1, ['mean', 'std']]
Out[143]:
In [144]:
m, s = mu_err_table['mu_err'].ix[1, ['mean', 'std']]
In [145]:
X = np.random.normal(m, s, size=numObjectsPerBin.ix[1] )
print(X)
In [146]:
XVals = []
for i in range(1,len(mu_err_table)):
m, s = mu_err_table['mu_err'].ix[1, ['mean', 'std']]
# We will convert the numpy array to list, but that may not be necessary
X = np.random.normal(m, s, size=numObjectsPerBin.ix[i]).tolist()
XVals.append(X)
In [217]:
sum(map(len, XVals))
Out[217]:
In [158]:
x=np.array(list(map(len, XVals)))
y=np.float(sum(list(map(len, XVals))))
In [160]:
x/y
Out[160]:
In [191]:
z[:numObjectsPerBin.values[0]]
for i in range(0,len(numObjectsPerBin)):
z[:numObjectsPerBin.values[i]]+= .1
In [ ]:
In [185]:
plt.hist(z, bins=np.arange(0,1.2,.1))
Out[185]:
In [199]:
from astropy.cosmology import Planck15
In [203]:
Planck15.distmod(z).value
Out[203]:
In [213]:
mu_err_or_something=np.concatenate(XVals)
In [214]:
len(mu_err_or_something)
Out[214]:
In [ ]:
In [ ]:
In [186]:
z
Out[186]:
In [190]:
z= np.random.uniform(0,.1, NumSN)
len(z)
Out[190]:
In [167]:
plt.hist(z)
Out[167]:
In [ ]:
In [31]:
mu_err_table.mu_err['count']
Out[31]:
In [32]:
def numinbins(tab,NumSN):
x=(tab.mu_err['count']/tab.mu_err['count'].sum())
return round(x*NumSN).astype(np.int)
In [33]:
numinbins(mu_err_table,2000)
Out[33]:
In [34]:
Out[34]:
In [ ]:
help (np.random.normal)
In [ ]:
help (np.histogram)
In [ ]:
hist, edges= np.histogram(z, bins=10)
In [ ]:
edges
In [ ]:
hist
In [ ]:
print (len(edges))
print (len(hist))
In [ ]:
binz=[]
for i in range(0,10):
stuff = edges[i]+edges[i+1]
point = stuff/2
binz.append (point)
print (binz)
In [ ]:
plt.plot(binz,hist)
In [ ]:
help (np.random.choice)
In [ ]:
wut=[]
for i in range(0,10):
k = hist[i]/hist.sum()
wut.append (k)
print (wut)
In [ ]:
print(len(binz))
print(len(wut))
In [ ]:
supernovastuff = np.random.choice(binz,100,p=wut)
print (supernovastuff)
In [ ]:
hist, edges= np.histogram(supernovastuff, bins=10)
In [ ]:
binzz=[]
for i in range(0,10):
stufff = edges[i]+edges[i+1]
pointt = stufff/2
binzz.append (pointt)
print (binzz)
In [ ]:
plt.plot(binzz,hist)
In [ ]:
supernovastuff = np.random.choice(binz,100,p=wut)
print (supernovastuff)
hist, edges= np.histogram(supernovastuff, bins=10)
binzz=[]
for i in range(0,10):
stufff = edges[i]+edges[i+1]
pointt = stufff/2
binzz.append (pointt)
plt.plot(binzz,hist)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: