In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from collections import OrderedDict as odict
import pandas as pd
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 [2]:
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[2]:
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 [3]:
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 [4]:
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 [5]:
#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 [6]:
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 [7]:
const = np.median(muz-mu)
print(const)
plt.scatter(z,muz-mu)
plt.plot([0,1.4],[const,const],'r')
Out[7]:
In [8]:
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 [9]:
help (np.histogram)
In [10]:
hist , edges = np.histogram(z, bins =10)
In [11]:
edges
Out[11]:
In [12]:
hist
Out[12]:
In [13]:
print(len (edges))
print(len (hist))
In [14]:
mid_pts =[]
for ii in range(len(edges)-1):
mid_pts.append(0.5*(edges[ii]+edges[ii+1]))
print(mid_pts)
In [ ]:
In [15]:
plt.plot (mid_pts, hist)
Out[15]:
In [16]:
help(np.random.choice)
np.random.choice(binz, 1000, p= #generate a random set of supernovae and see if we an get a different set of random z values and see if we can get them to look the same) # number of supernova in the bin vs total number of supernova
In [ ]:
In [17]:
plt.plot(z, np.sqrt(mu_var), 'o')
Out[17]:
In [18]:
import pandas as pd
In [19]:
df = pd.dataFrame(gzip.GzipFile('snFits.p.gz'),
encoding='latin1')
In [ ]:
type (snFits)
In [20]:
df = pd.DataFrame(snFits)
In [21]:
df.head()
Out[21]:
In [22]:
snf = [s for s in snFits.values() if s is not None]
In [23]:
len(snf)
Out[23]:
In [24]:
snFits['18263']
Out[24]:
In [25]:
df = df.transpose()
In [26]:
results = df[['mu', 'mu_var', 'z']].astype(np.float)
results.head()
Out[26]:
In [27]:
results.dtypes
Out[27]:
In [28]:
lowz = results.query('z < .4')
lowz.head()
Out[28]:
In [29]:
plt.hist(lowz['mu_var'], normed=1)
Out[29]:
In [30]:
plt.hist (highz['mu_var'],normed=1 )
In [31]:
results['bindex'] = results['z'] // 0.1
In [32]:
results.head()
Out[32]:
In [33]:
results['mu_err'] = np.sqrt(results['mu_var'])
results.head()
Out[33]:
In [34]:
results.dtypes
Out[34]:
In [ ]:
In [35]:
grouped = results.groupby('bindex')
In [36]:
help(.agg)
In [ ]:
mu_error= np.sqrt(mu_var)
In [37]:
def clipped_mean(mu_err, outlier_threshhold):
return(mu_err[mu_err< outlier_threshhold])
In [38]:
mu_err_table = grouped.agg(dict(mu_err=[np.mean, np.std, "count"]))
mu_err_table
Out[38]:
In [39]:
mu_err_table.mu_err['count'].sum()
Out[39]:
In [40]:
def numinbins(tab, NumSN):
x = (tab.mu_err['count']/tab.mu_err['count'].sum())
return(round(x*NumSN).astype(np.int))
In [51]:
numinbins(mu_err_table, 2000)
Out[51]:
In [42]:
#write muerr
#loop through each bin and take the number we get from num in bins and then
#np.random.normal(mean,std, size)
In [47]:
def muerr(tab, NumSN):
return np.random.normal(tab.mu_err['mean'], tab.mu_err ['std'])
muerr(mu_err_table, 1000)
Out[47]:
In [ ]:
In [ ]:
In [ ]:
In [52]:
mu_err_table['frequencies']= mu_err_table.mu_err['count']/mu_err_table.mu_err['count'].sum()
In [69]:
mu_err_table= mu_err_table.ix[1:]
In [70]:
mu_err_table
Out[70]:
mu_err_table
In [100]:
NumSN= 50000
numObjectsPerBin = np.round(mu_err_table.frequencies * NumSN).astype(np.int)
print(numObjectsPerBin)
In [ ]:
In [78]:
m, s = [mu_err_table.mu_err['mean'], mu_err_table.mu_err['std']] # Now the mean of the 0th bin is assigned to m, and std to s
#X = np.random.normal(m, s, size=numObjectsPerBin.ix[1] )
#print()
In [79]:
m
Out[79]:
In [80]:
s
Out[80]:
In [87]:
numObjectsPerBin
Out[87]:
In [134]:
X = [np.random.normal(m[i], s[i], size=numObjectsPerBin[i]) for i in range(1, len(numObjectsPerBin))]
In [135]:
X[0]
Out[135]:
In [90]:
len(X[0])
Out[90]:
In [91]:
np.mean(X[0])
Out[91]:
In [92]:
np.mean(X[3])
Out[92]:
In [133]:
len(X)#x=xvals
Out[133]:
In [94]:
type(X)
Out[94]:
In [97]:
x=np.array(list(map(len, X)))
y =np.float(sum(list(map(len,X))))
In [103]:
z= np.random.uniform(0,.1, size = NumSN)
len(z)
Out[103]:
In [104]:
plt.hist(z)
Out[104]:
In [105]:
numObjectsPerBin
Out[105]:
In [112]:
z[:numObjectsPerBin.values[0]]
Out[112]:
In [113]:
z[:numObjectsPerBin.values[i]]
Out[113]:
In [116]:
for i in range (len(numObjectsPerBin)):
z[:numObjectsPerBin.values[i]]+= .1
In [121]:
plt.hist(z, bins= np.arange(0,1.2,.1))
Out[121]:
In [122]:
z[10:]
Out[122]:
In [123]:
z
Out[123]:
In [127]:
from astropy.cosmology import Planck15
In [128]:
Planck15.distmod(z).value
Out[128]:
In [136]:
X
Out[136]:
In [ ]: