In [5]:
%matplotlib inline
import matplotlib.pyplot as plt
import os
import numpy as np
from astropy.units import Unit
from astropy.table import Table, join
import collections
import pandas as pd
import sncosmo
import triangle
import lc
In [6]:
bandPassList = ['u', 'g', 'r', 'i', 'z', 'y']
banddir = os.path.join(os.getenv('THROUGHPUTS_DIR'), 'baseline')
# lsstbands = list()
# lsstbp = dict()
for band in bandPassList:
# setup sncosmo bandpasses
bandfname = banddir + "/total_" + band + '.dat'
# register the LSST bands to the SNCosmo registry
# Not needed for LSST, but useful to compare independent codes
# Usually the next two lines can be merged,
# but there is an astropy bug currently which affects only OSX.
numpyband = np.loadtxt(bandfname)
sncosmoband = sncosmo.Bandpass(wave=numpyband[:, 0],
trans=numpyband[:, 1],
wave_unit=Unit('nm'),
name='LSST_' + band)
sncosmo.registry.register(sncosmoband, force=True)
In [9]:
model = sncosmo.Model(source='salt2-extended')
params = {'z': 1.06371, 't0': 50916.5, 'x0': 1.695e-06, 'x1': -0.708466, 'c': 0.0178018}
model.set(**params)
vparams = ['t0', 'x0', 'x1', 'c']
In [11]:
test = lc.LC(model, vparams)
In [12]:
test.fitOut
Out[12]:
In [13]:
test.mcmcOut
Out[13]:
In [8]:
test.nestOut
Out[8]:
In [231]:
test.mcmcRes.samples.shape
Out[231]:
In [28]:
# syntactic sugar
fit_errors = test.fitRes.errors
In [29]:
fit_table = Table()
for keys in fit_errors:
fit_table[keys] = [fit_errors[keys]]
for key in test.fitRes.keys():
if key == 'errors':
continue
fit_table.meta[key] = test.fitRes[key]
In [30]:
fit_table, fit_table.meta
Out[30]:
In [173]:
fit_table.write('Fits.hdf5', 'MLEfit')
In [174]:
# syntactic sugar
mcmc_errors = test.mcmcRes.errors
mcmc_table = Table(test.mcmcRes.samples, names=test.mcmcRes.vparam_names)
for key in test.mcmcRes.keys():
if key == 'errors' or key =='samples':
continue
mcmc_table.meta[key] = test.mcmcRes[key]
mcmc_table.add_row(mcmc_errors.values())
mcmc_table, mcmc_table.meta
Out[174]:
In [175]:
mcmc_table.write('Fits.hdf5', 'mcmc', append=True)
In [176]:
test.nestRes
Out[176]:
In [177]:
# syntactic sugar
nest_errors = test.nestRes.errors
nest_param_dict = test.nestRes.param_dict
nest_bounds = test.nestRes.bounds
nest_table = Table(test.nestRes.samples, names=test.nestRes.vparam_names)
In [179]:
for key in test.nestRes.keys():
if key == 'errors' or key =='samples' or key =='param_dict' or key == 'bounds':
continue
nest_table.meta[key] = test.nestRes[key]
nest_table.add_row(nest_errors.values())
temp_z = np.zeros((len(nest_table['t0'])))
col_z = Table.Column(name='z', data=temp_z)
nest_table.add_column(col_z)
In [180]:
param_list = []
for colname in nest_table.colnames:
param_list.append(nest_param_dict[colname])
nest_table.add_row(param_list)
In [182]:
for key in nest_bounds.keys():
nest_table.meta[key] = nest_bounds[key]
In [183]:
nest_table.write('Fits.hdf5', 'nest', append=True)
In [184]:
fit_read = Table.read('Fits.hdf5', 'MLEfit')
fit_read, fit_read.meta
Out[184]:
In [185]:
fit_errors_dict = collections.OrderedDict()
for colnames in fit_read.colnames:
fit_errors_dict[colnames] = fit_read[colnames][0]
In [186]:
fit_dict = fit_read.meta
fit_dict['errors'] = fit_errors_dict
Out[186]:
In [187]:
fit_result = sncosmo.utils.Result(fit_dict)
fit_result
Out[187]:
In [188]:
mcmc_read = Table.read('Fits.hdf5', 'mcmc')
Out[188]:
In [189]:
mcmc_errors_dict = collections.OrderedDict()
for colnames in mcmc_read.colnames:
mcmc_errors_dict[colnames] = mcmc_read[colnames][len(mcmc_read.columns[0]) - 1]
In [190]:
mcmc_read.remove_row(len(mcmc_read.columns[0]) - 1)
In [255]:
mcmc_samples = np.array([np.array(mcmc_read.columns[0]),
np.array(mcmc_read.columns[1]),
np.array(mcmc_read.columns[2]),
np.array(mcmc_read.columns[3])])
In [256]:
mcmc_dict = mcmc_read.meta
mcmc_dict['errors'] = mcmc_errors_dict
mcmc_dict['samples'] = mcmc_samples.T
mcmc_result = sncosmo.utils.Result(mcmc_dict)
mcmc_result
Out[256]:
In [257]:
nest_read = Table.read('Fits.hdf5', 'nest')
In [258]:
nest_param_dict = collections.OrderedDict()
for colnames in nest_read.colnames:
nest_param_dict[colnames] = nest_read[colnames][len(nest_read.columns[0]) - 1]
nest_read.remove_row(len(nest_read.columns[0]) - 1)
nest_read.remove_column('z')
nest_errors_dict = collections.OrderedDict()
for colnames in nest_read.colnames:
nest_errors_dict[colnames] = nest_read[colnames][len(nest_read.columns[0]) - 1]
nest_read.remove_row(len(nest_read.columns[0]) - 1)
In [259]:
nest_samples = np.array([np.array(nest_read.columns[0]),
np.array(nest_read.columns[1]),
np.array(nest_read.columns[2]),
np.array(nest_read.columns[3])])
In [260]:
nest_bounds = {}
for colnames in nest_read.colnames:
nest_bounds[colnames] = tuple(nest_read.meta[colnames])
del nest_read.meta[colnames]
In [261]:
nest_dict = nest_read.meta
nest_dict['errors'] = nest_errors_dict
nest_dict['param_dict'] = nest_param_dict
nest_dict['samples'] = nest_samples.T
nest_dict['bounds'] = nest_bounds
In [262]:
nest_result = sncosmo.utils.Result(nest_dict)
nest_result
Out[262]:
In [263]:
io_check = lc.LC(model)
In [264]:
io_check.fitRes = fit_result
io_check.mcmcRes = mcmc_result
io_check.nestRes = nest_result
In [265]:
io_check.fitRes
Out[265]:
In [269]:
fit_model_dict = {}
index = 0
for vals in io_check.fitRes.param_names:
fit_model_dict[vals] = io_check.fitRes.parameters[index]
index = index + 1
In [270]:
mcmc_model_dict = {}
index = 0
for vals in io_check.mcmcRes.param_names:
mcmc_model_dict[vals] = io_check.mcmcRes.parameters[index]
index = index + 1
In [271]:
nest_model_dict = io_check.nestRes.param_dict
In [272]:
fitmodel = sncosmo.Model(source='salt2-extended')
fitmodel.set(**fit_model_dict)
io_check.fitModel = fitmodel
print io_check.fitModel
In [273]:
mcmcmodel = sncosmo.Model(source='salt2-extended')
mcmcmodel.set(**mcmc_model_dict)
io_check.mcmcModel = mcmcmodel
print io_check.mcmcModel
In [274]:
nestmodel = sncosmo.Model(source='salt2-extended')
nestmodel.set(**nest_model_dict)
io_check.nestModel = nestmodel
print io_check.nestModel
In [275]:
fig = lc.LC.plotLC(io_check)
In [276]:
lc.LC.plotCorner(io_check)
Out[276]:
In [277]:
lc.LC.plotTrace(io_check)
Out[277]:
In [ ]: