In [1]:
from LSSTmetrics import PerSNMetric, EfficiencyTable
from LSSTmetrics import example_data as lsstm_example_data
from opsimsummary import OpSimOutput
#from efficiencyTable import EfficiencyTable
In [2]:
from opsimsummary import summarize_opsim as oss
In [3]:
import pandas as pd
import sncosmo
In [4]:
from sndata import LightCurve
In [5]:
%matplotlib inline
import matplotlib.pyplot as plt
import os
import numpy as np
import copy
In [6]:
from lsst.sims.photUtils import BandpassDict
In [7]:
# Catsim bandpasses
lsst_bp = BandpassDict.loadTotalBandpassesFromFiles()
In [8]:
# sncosmo Bandpasses required for fitting
throughputsdir = os.getenv('THROUGHPUTS_DIR')
from astropy.units import Unit
bandPassList = ['u', 'g', 'r', 'i', 'z', 'y']
banddir = os.path.join(os.getenv('THROUGHPUTS_DIR'), 'baseline')
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)
print band
sncosmoband = sncosmo.Bandpass(wave=numpyband[:, 0],
trans=numpyband[:, 1],
wave_unit=Unit('nm'),
name=band)
sncosmo.registry.register(sncosmoband, force=True)
In [9]:
efficiencyFile = os.path.join(lsstm_example_data, 'SEARCHEFF_PIPELINE_DES.DAT')
et = EfficiencyTable.fromDES_EfficiencyFile(efficiencyFile)
In [10]:
from lsst.sims.catUtils.supernovae import SNObject
In [11]:
opsimHDF = os.path.join(os.getenv('HOME'), 'data', 'LSST', 'OpSimData', 'storage.h5')
#summarydf = OpSimOutput.fromOpSimDB('/Users/rbiswas/data/LSST/OpSimData/minion_1016_sqlite.db').summary.copy()
summarydf = pd.read_hdf(opsimHDF, 'table')
# df = df.query('propID == [364, 366]')
In [12]:
summarydf.head()
Out[12]:
In [13]:
#summarydf = summarydf.query('propID == [54, 56]')#.query('night < 365')
In [14]:
# Create the summary instance
so = oss.SummaryOpsim(summarydf)
In [15]:
ss = so.simlib(fieldID=309)
In [16]:
ss.head()
Out[16]:
In [17]:
fig_309, mat, _, _, _ = so.cadence_plot(fieldID=309, nightMax=730, nightMin=365., sql_query='night < 730. and night > 365.')
In [18]:
fig_744, mat, _, _, _ = so.cadence_plot(fieldID=744, nightMax=730, nightMin=365., sql_query='night < 730. and night > 365.')
In [19]:
fig_309.savefig('fig_309_2ndYear.pdf')
fig_744.savefig('fig_744_2ndYear.pdf')
In [20]:
so.mjdvalfornight(570)
Out[20]:
In [21]:
fig_744_49923, _, _, _, _ = so.cadence_plot(fieldID=744, mjd_center=49923,
mjd_range=[-20., 50.])
In [66]:
fig_309_49923, _, _, _, _ = so.cadence_plot(fieldID=309, mjd_center=49923,
mjd_range=[-20., 50.])
In [22]:
fig_309_49823, _, _, _, _ = so.cadence_plot(fieldID=309, mjd_center=49823,
mjd_range=[-20., 50.])
In [24]:
fig_309_49540, _, _, _, _ = so.cadence_plot(fieldID=309, mjd_center=49540,
mjd_range=[-20., 50.])
In [25]:
fig_290_49540, _, _, _, _ = so.cadence_plot(fieldID=290, mjd_center=49540,
mjd_range=[-20., 50.])
In [30]:
fig_290_49540.savefig('SN_Cadence_290.pdf')
In [29]:
fig_309_49540.savefig('SN_Cadence_309.pdf')
In [79]:
fig_744_49923.savefig('TimeWindow_744_49923.pdf')
In [75]:
fig_309_49923.savefig('TimeWindow_309_49923.pdf')
In [76]:
fig_309_49823.savefig('TimeWindow_309_49823.pdf')
In [31]:
fig_309_49540, _, _, _, _ = so.cadence_plot(fieldID=309, mjd_center=49540,
mjd_range=[-20., 50.])
In [32]:
qm_290_49540 = PerSNMetric(t0=49540, fieldID=290, summarydf=summarydf, efficiency=et,
lsst_bp=lsst_bp)
In [33]:
qm_290_49540.lightcurve
Out[33]:
In [34]:
import seaborn as sns
sns.set_style('whitegrid')
sns.set_context('paper')
In [35]:
qm_309_49540 = PerSNMetric(t0=49540, fieldID=309, summarydf=summarydf, efficiency=et,
lsst_bp=lsst_bp)
In [36]:
def plotlc(qm, ylabel='flux(maggies)', xlabel_add='(MJD)'):
model = qm.sncosmoModel
lc = LightCurve(qm.lightcurve)
fig = sncosmo.plot_lc(lc.snCosmoLC(), model=model, color='k', pulls=False)
axs = fig.axes
xlabel = axs[-1].get_xlabel()
xlabel += xlabel_add
_ = list(ax.set_ylabel(ylabel) for ax in axs[::2])
_ = list(ax.set_xlabel(xlabel) for ax in axs[-2:])
return fig
In [37]:
sns.set_style(None)
In [38]:
fig_309_lc = plotlc(qm_309_49540)
In [39]:
fig_309_lc.savefig('SN_309_lc.pdf')
In [40]:
fig_290_lc = plotlc(qm_290_49540)
In [41]:
fig_290_lc.axes[2].get_ylabel()
Out[41]:
In [42]:
fig_290_lc.savefig('SN_290_lc.pdf')
In [43]:
qm_744_499 = PerSNMetric(t0=49923, fieldID=744, summarydf=summarydf, efficiency=et,
lsst_bp=lsst_bp)
In [44]:
qm_309_499 = PerSNMetric(t0=49923, fieldID=309, summarydf=summarydf, efficiency=et,
lsst_bp=lsst_bp)
In [45]:
qm_309_498 = PerSNMetric(t0=49823, fieldID=309, summarydf=summarydf, efficiency=et,
lsst_bp=lsst_bp)
In [59]:
fig_744_499_lc = qm_744_499.lcplot()
fig_309_499_lc = qm_309_499.lcplot()
fig_309_498_lc = qm_309_498.lcplot()
# save to files
fig_744_499_lc.savefig('fig_744_49923_lc.pdf')
fig_309_499_lc.savefig('fig_309_49923_lc.pdf')
fig_309_498_lc.savefig('fig_309_49823_lc.pdf')
In [60]:
mcmc_out_309_744_Nocov = sncosmo.mcmc_lc(qm_744_499.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=False)
mcmc_out_309_499_Nocov = sncosmo.mcmc_lc(qm_309_499.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=False)
mcmc_out_309_498_Nocov = sncosmo.mcmc_lc(qm_309_498.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=False)
In [ ]:
qm_309_499.writeLightCurve('lc_field309_mjd_49923.ascii')
In [ ]:
qm_744_499.writeLightCurve('lc_field744_mjd_49923.ascii')
In [24]:
fig_744_499_lc
Out[24]:
In [214]:
res_309_498_NoCov = ans.ResChar.fromSNCosmoRes(mcmc_out_309_498_Nocov)
In [218]:
corner(res_309_498_NoCov.salt_samples(), labels=res_309_498_NoCov.salt_samples().columns,
truths=[49823., 0., 0., -12.5, -12.5])
Out[218]:
In [ ]:
res_744_499_NoCov = ans.ResCharrChar.fromSN
In [ ]:
mcmc_out_309_498 = sncosmo.mcmc_lc(qm_309_498.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=True)
In [193]:
qm_309_498 = PerSNMetric(t0=49823, fieldID=309, summarydf=summarydf, efficiency=et,
lsst_bp=lsst_bp)
In [198]:
fig_qm_309_498 = qm_309_498.lcplot()
In [206]:
fig_744_499 = qm_744_499.lcplot(nightlyCoadd=False)
In [ ]:
In [183]:
qm_309_498.SN.SNstate
Out[183]:
In [185]:
print qm_309_498.SN.equivalentSNCosmoModel()
In [156]:
mcmc_out_309_498 = sncosmo.mcmc_lc(qm_309_498.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=True)
In [157]:
out_309_498 = ans.ResChar.fromSNCosmoRes(mcmc_out_309_498)
In [162]:
-2.5 * np.log10(qm_309_498.SN.get('x0'))
Out[162]:
In [164]:
qm_309_498.lcplot()
Out[164]:
In [181]:
print qm_309_498.SN.equivalentSNCosmoModel()
In [182]:
print out_309_498.sncosmoModel
In [186]:
sncosmo.chisq(qm_309_498.SNCosmoLC(), model=qm_309_498.SN.equivalentSNCosmoModel())
Out[186]:
In [191]:
sncosmo.chisq(qm_309_498.SNCosmoLC(), model=out_309_498.sncosmoModel)
Out[191]:
In [190]:
print qm_309_498.sncosmoModel
In [176]:
sncosmo.plot_lc( qm_309_498.SNCosmoLC(), model=[out_309_498.sncosmoModel,
qm_309_498.SN.equivalentSNCosmoModel()],
model_label=['fitted', 'truth'])
Out[176]:
In [163]:
corner(out_309_498.salt_samples(), labels=out_309_498.salt_samples().columns, truths=[49823, 12.5, 0, 0., 12.5])
Out[163]:
In [80]:
qm_744 = PerSNMetric(t0=49923, fieldID=744, summarydf=summarydf, efficiency=et, lsst_bp=lsst_bp)
In [82]:
qm_744.qualityMetric()
Out[82]:
In [139]:
mBT = -2.5 * np.log10(qm_744.SN.SNstate['x0'])
muT = mBT
In [88]:
import analyzeSN as ans
from corner import corner
In [126]:
mcmc_out_744_F = sncosmo.mcmc_lc(qm_744.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=False)
In [127]:
tt = ans.ResChar.fromSNCosmoRes(mcmc_out_744_F)
In [129]:
tt.salt_samples().mu.var()/ out_744.salt_samples().mu.var()
Out[129]:
In [90]:
mcmc_out = sncosmo.mcmc_lc(qm_744.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=True)
In [120]:
mcmc_out_309 = sncosmo.mcmc_lc(qm_309.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=True)
In [130]:
mcmc_out_309_F = sncosmo.mcmc_lc(qm_309.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=False)
In [131]:
out_309_F = ans.ResChar.fromSNCosmoRes(mcmc_out_309_F)
In [91]:
out_744 = ans.ResChar.fromSNCosmoRes(mcmc_out)
In [122]:
out_309 = ans.ResChar.fromSNCosmoRes(mcmc_out_309)
In [132]:
out_309.salt_samples().mu.var() / out_309_F.salt_samples().mu.var()
Out[132]:
In [124]:
0.05 * 0.05 / out_309.salt_samples().mu.var()
Out[124]:
In [125]:
0.05 * 0.05 / out_309.mu_variance_linear()
Out[125]:
In [142]:
qm_309.lcplot(scattered=False)
Out[142]:
In [143]:
corner(out_309.salt_samples(), labels=out_309.salt_samples().columns, truths=[49923,
0., 0. , mBT, muT ])
Out[143]:
In [94]:
corner(out_744.salt_samples(), labels=out_744.salt_samples().columns)
Out[94]:
In [115]:
0.05**2 / out_744.salt_samples(alpha=0.11, beta=3.1).mu.var()
Out[115]:
In [119]:
0.05**2 / out_fit_309.salt_samples(alpha=0.11, beta=3.1).mu.var()
In [113]:
out_fit_744 = sncosmo.fit_lc(qm_744.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['x0', 'x1','c'], minsnr=0.,)
In [117]:
out_fit_744 = sncosmo.fit_lc(qm_744.SNCosmoLC(nightlyCoadd=True), model=qm_744.sncosmoModel,
vparam_names=['x0', 'x1','c'], minsnr=0.,)
In [114]:
ans.ResChar.fromSNCosmoRes(out_fit_744)
In [105]:
qm_744.discoveryMetric()
Out[105]:
In [78]:
qm_309 = PerSNMetric(t0=49923, fieldID=309, summarydf=summarydf, efficiency=et, lsst_bp=lsst_bp)
In [79]:
qm_309.discoveryMetric()
Out[79]:
In [84]:
qm_309.qualityMetric()
Out[84]:
In [51]:
loc, ax_loc =plt.subplots()
In [57]:
xx = oss.SummaryOpsim(summarydf.query('fieldID==309'))
loc_fig = xx.showFields()
In [58]:
yy = oss.SummaryOpsim(summarydf.query('fieldID==744'))
yy.showFields(ax=loc_fig.axes[0])
Out[58]:
In [60]:
loc_fig.savefig('loc_309_744.pdf')
In [37]:
mat[mat.fillna(0.).sum(axis=1) > 0].index.max() - mat[mat.fillna(0.).sum(axis=1) > 0].index.min()
Out[37]:
In [24]:
len(_)
Out[24]:
In [16]:
_ = so.cadence_plot(fieldID=309, mjd_center=49540)
In [17]:
# 1st season by default
_ = so.cadence_plot(fieldID=309)
In [80]:
qm = PerSNMetric(fieldID=309, t0=49540, summarydf=summarydf, lsst_bp=lsst_bp)
In [82]:
fig_49540_no_coadd = qm.lcplot()
In [83]:
fig_49540_no_coadd.axes
Out[83]:
In [63]:
qm.lcplot(nightlyCoadd=True)
Out[63]:
In [21]:
qm.writeLightCurve('coadded_lc.dat',nightlyCoadd=True)
In [22]:
cm = qm.cadence_Matrix(fieldID=309, mjd_center=49540, mjd_range=[-30., 50.])
In [23]:
cm.head()
Out[23]:
In [24]:
cm = qm.cadence_plot(fieldID=309, mjd_center=49540, mjd_range=[-30, 5])
In [25]:
qm = PerSNMetric(fieldID=309, t0=49540, summarydf=ss, lsst_bp=lsst_bp, efficiency=et)
In [26]:
qm.SNCadence[0]
Out[26]:
In [27]:
qm.summary.head()
Out[27]:
In [28]:
qm.discoveryMetric()
Out[28]:
In [29]:
qm.qualityMetric()
Out[29]:
In [30]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import healpy as hp
import lsst.sims.maf.db as db
import lsst.sims.maf.utils as utils
import lsst.sims.maf.metrics as metrics
import lsst.sims.maf.slicers as slicers
import lsst.sims.maf.metricBundles as metricBundles
In [31]:
outDir ='LightCurve'
dbFile = 'enigma_1189_sqlite.db'
opsimdb = utils.connectOpsimDb(dbFile)
resultsDb = db.ResultsDb(outDir=outDir)
In [32]:
filters = ['u','g','r','i','z','y']
colors={'u':'cyan','g':'g','r':'y','i':'r','z':'m', 'y':'k'}
In [33]:
# Set RA, Dec for a single point in the sky. in radians.
ra = np.radians(0.)
dec = np.radians(0.)
# SNR limit (Don't use points below this limit)
snrLimit = 5.
# Demand this many points above SNR limit before plotting LC
nPtsLimit = 6
In [34]:
# The pass metric just passes data straight through.
metric = metrics.PassMetric(cols=['filter','fieldID','finSeeing','fiveSigmaDepth',
'expMJD','airmass', 'propID', 'night', 'filtSkyBrightness'])
slicer = slicers.UserPointsSlicer(ra,dec,lonCol='fieldRA',latCol='fieldDec')
sql = 'night < 366'
bundle = metricBundles.MetricBundle(metric,slicer,sql)
bg = metricBundles.MetricBundleGroup({0:bundle}, opsimdb,
outDir=outDir, resultsDb=resultsDb)
In [35]:
bg.runAll()
In [36]:
llc = pd.DataFrame.from_records(bundle.metricValues.data[0])
In [37]:
llc.head()
Out[37]:
In [38]:
llc.expMJD.hist()
Out[38]:
In [39]:
q2 = PerSNMetric(t0=49580, summarydf=llc, lsst_bp=lsst_bp, efficiency=et, raCol='fieldRA',
decCol='fieldDec')
In [40]:
ss.columns
Out[40]:
In [41]:
llc.columns
Out[41]:
In [42]:
q2.raCol
Out[42]:
In [43]:
q2.summary[q2.raCol].iloc[0]
Out[43]:
In [44]:
q2.radeg
Out[44]:
In [45]:
q2.decdeg
Out[45]:
In [46]:
q2.SN.SNstate
Out[46]:
In [47]:
q2.lcplot(nightlyCoadd=False)
In [48]:
q2.lcplot(nightlyCoadd=True)
In [49]:
q2.writeLightCurve(fname='test.dat', nightlyCoadd=True)
In [50]:
!cat test.dat
In [51]:
q2.coaddedLightCurve.reset_index()[q2.coaddedLightCurve.reset_index()['band'] == 'y']
Out[51]:
In [52]:
q2.qualityMetric()
Out[52]:
In [53]:
q2.discoveryMetric()
Out[53]:
In [54]:
_ = q2.cadence_plot(summarydf=llc, racol='fieldRA', deccol='fieldDec', mjd_center=49580)
In [55]:
# Defaults to full season
_ = q2.cadence_plot(summarydf=llc, racol='fieldRA', deccol='fieldDec')
In [56]:
xx = llc[['expMJD', 'filter', 'fiveSigmaDepth', 'fieldRA', 'fieldDec', 'fieldID']].copy(deep=True)
In [57]:
xx
Out[57]:
In [58]:
q3 = PerSNMetric(t0=49580, summarydf=xx, lsst_bp=lsst_bp, efficiency=et, raCol='fieldRA',
decCol='fieldDec')
In [59]:
q3.lightcurve
Out[59]:
In [60]:
q3.SNCosmoLC()
Out[60]:
In [61]:
q3.SN.SNstate
Out[61]:
In [62]:
q3.qualityMetric()
Out[62]:
In [63]:
q3.discoveryMetric()
Out[63]:
In [ ]: