In [1]:
import numpy as np

In [2]:
import sncosmo

In [3]:
%matplotlib inline
import matplotlib.pyplot as plt

In [ ]:


In [4]:
import sncosmo

In [5]:
from analyzeSN import LightCurve

In [6]:
lco = LightCurve.fromSALTFormat('/Users/rbiswas/data/JLA_Calibration_OnlineData/cosmology/jla_light_curves/lc-03D1au.list')

In [7]:
sncosmo.plot_lc(lco.snCosmoLC())


Out[7]:

In [8]:
lco = LightCurve.fromSALTFormat('/Users/rbiswas/data/JLA_Calibration_OnlineData/cosmology/jla_light_curves/lc-SDSS14331.list')

In [9]:
sncosmo.plot_lc(lco.snCosmoLC())


Out[9]:

In [9]:
t = sncosmo.read_lc('/Users/rbiswas/data/JLA_Calibration_OnlineData/cosmology/jla_light_curves/lc-SDSS14331.list',
                   format='salt2')

In [10]:
t


Out[10]:
<Table length=106>
DateFluxFluxerrZPFilterMagSys
float64float64float64float64str7str6
53974.332-3.6712.33727.5SDSS::gAB_B12
53974.332-3.7718.60827.5SDSS::rAB_B12
53974.332-13.3929.72527.5SDSS::iAB_B12
53974.3325.9991.96327.5SDSS::zAB_B12
53974.3326.9737.97827.5SDSS::uAB_B12
53977.383-13.9612.59927.5SDSS::gAB_B12
53977.383-0.619.24827.5SDSS::rAB_B12
53977.38335.2530.52627.5SDSS::iAB_B12
53977.38317.0796.46227.5SDSS::zAB_B12
53977.38321.4636.40427.5SDSS::uAB_B12
..................
54060.2811.2314.61427.5SDSS::gAB_B12
54060.28169.127.79327.5SDSS::rAB_B12
54060.28119.1344.18827.5SDSS::iAB_B12
54060.28163.99129.7227.5SDSS::zAB_B12
54060.281-8.0953.76527.5SDSS::uAB_B12
54063.2517.6613.94927.5SDSS::gAB_B12
54063.2535.8722.08627.5SDSS::rAB_B12
54063.2579.635.71427.5SDSS::iAB_B12
54063.25141.14122.87527.5SDSS::zAB_B12
54063.251.1544.68427.5SDSS::uAB_B12

In [8]:
lco._lightCurve


Out[8]:
Date Flux Fluxerr ZP Filter MagSys
0 53974.332 -3.67 12.337 27.5 SDSS::g AB_B12
1 53974.332 -3.77 18.608 27.5 SDSS::r AB_B12
2 53974.332 -13.39 29.725 27.5 SDSS::i AB_B12
3 53974.332 5.99 91.963 27.5 SDSS::z AB_B12
4 53974.332 6.97 37.978 27.5 SDSS::u AB_B12
5 53977.383 -13.96 12.599 27.5 SDSS::g AB_B12
6 53977.383 -0.60 19.248 27.5 SDSS::r AB_B12
7 53977.383 35.25 30.526 27.5 SDSS::i AB_B12
8 53977.383 17.07 96.462 27.5 SDSS::z AB_B12
9 53977.383 21.46 36.404 27.5 SDSS::u AB_B12
10 53989.336 15.85 30.719 27.5 SDSS::g AB_B12
11 53989.336 -5.17 35.300 27.5 SDSS::r AB_B12
12 53989.336 3.05 39.413 27.5 SDSS::i AB_B12
13 53989.336 -3.01 100.574 27.5 SDSS::z AB_B12
14 53997.434 69.58 16.132 27.5 SDSS::g AB_B12
15 53997.434 71.77 21.879 27.5 SDSS::r AB_B12
16 54000.313 116.15 17.035 27.5 SDSS::g AB_B12
17 54000.313 123.24 24.404 27.5 SDSS::r AB_B12
18 54000.313 140.42 38.033 27.5 SDSS::i AB_B12
19 54000.313 165.47 141.974 27.5 SDSS::z AB_B12
20 54000.313 -15.25 49.873 27.5 SDSS::u AB_B12
21 54006.258 313.65 20.874 27.5 SDSS::g AB_B12
22 54006.258 290.85 26.311 27.5 SDSS::r AB_B12
23 54006.258 1590.24 257.563 27.5 SDSS::z AB_B12
24 54007.363 291.60 16.527 27.5 SDSS::g AB_B12
25 54007.363 318.41 21.851 27.5 SDSS::r AB_B12
26 54007.363 337.89 33.644 27.5 SDSS::i AB_B12
27 54007.363 269.70 112.001 27.5 SDSS::z AB_B12
28 54007.363 49.11 43.277 27.5 SDSS::u AB_B12
29 54009.285 337.36 17.824 27.5 SDSS::g AB_B12
... ... ... ... ... ... ...
76 54052.266 40.84 24.873 27.5 SDSS::g AB_B12
77 54052.266 38.61 37.564 27.5 SDSS::r AB_B12
78 54052.266 33.42 58.622 27.5 SDSS::i AB_B12
79 54052.266 105.99 205.316 27.5 SDSS::z AB_B12
80 54052.266 88.97 64.142 27.5 SDSS::u AB_B12
81 54054.285 6.61 25.251 27.5 SDSS::g AB_B12
82 54054.285 64.02 41.704 27.5 SDSS::r AB_B12
83 54054.285 245.29 56.497 27.5 SDSS::i AB_B12
84 54054.285 101.99 200.486 27.5 SDSS::z AB_B12
85 54054.285 159.18 84.760 27.5 SDSS::u AB_B12
86 54055.262 20.19 15.100 27.5 SDSS::g AB_B12
87 54055.262 59.86 22.113 27.5 SDSS::r AB_B12
88 54055.262 82.21 32.209 27.5 SDSS::i AB_B12
89 54055.262 69.37 106.812 27.5 SDSS::z AB_B12
90 54055.262 55.73 47.638 27.5 SDSS::u AB_B12
91 54058.262 74.93 16.635 27.5 SDSS::g AB_B12
92 54058.262 39.11 26.893 27.5 SDSS::r AB_B12
93 54058.262 151.06 41.455 27.5 SDSS::i AB_B12
94 54058.262 146.65 146.087 27.5 SDSS::z AB_B12
95 54058.262 7.21 52.412 27.5 SDSS::u AB_B12
96 54060.281 1.23 14.614 27.5 SDSS::g AB_B12
97 54060.281 69.10 27.793 27.5 SDSS::r AB_B12
98 54060.281 19.13 44.188 27.5 SDSS::i AB_B12
99 54060.281 63.99 129.720 27.5 SDSS::z AB_B12
100 54060.281 -8.09 53.765 27.5 SDSS::u AB_B12
101 54063.250 17.66 13.949 27.5 SDSS::g AB_B12
102 54063.250 35.87 22.086 27.5 SDSS::r AB_B12
103 54063.250 79.60 35.714 27.5 SDSS::i AB_B12
104 54063.250 141.14 122.875 27.5 SDSS::z AB_B12
105 54063.250 1.15 44.684 27.5 SDSS::u AB_B12

106 rows × 6 columns


In [7]:
lco.lightCurve.head()


Out[7]:
mjd flux fluxerr zp band zpsys
0 53974.332 -3.67 12.337 27.5 sdssg AB_B12
1 53974.332 -3.77 18.608 27.5 sdssr AB_B12
2 53974.332 -13.39 29.725 27.5 sdssi AB_B12
3 53974.332 5.99 91.963 27.5 sdssz AB_B12
4 53974.332 6.97 37.978 27.5 sdssu AB_B12

In [5]:
sncosmo.plot_lc(lco.snCosmoLC())`


---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-5-0d08b6d33e8e> in <module>()
----> 1 sncosmo.plot_lc(lco.snCosmoLC())

/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/plotting.pyc in plot_lc(data, model, bands, zp, zpsys, pulls, xfigsize, yfigsize, figtext, model_label, errors, ncol, figtextsize, show_model_params, tighten_ylim, color, cmap, cmap_lims, fname, **kwargs)
    184     if data is not None:
    185         data = standardize_data(data)
--> 186         data = normalize_data(data, zp=zp, zpsys=zpsys)
    187 
    188     # Bands to plot

/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/photdata.pyc in normalize_data(data, zp, zpsys)
    154         for ms in set(bandzpsys):
    155             idx2 = bandzpsys == ms
--> 156             ms = get_magsystem(ms)
    157             bandfactor[idx2] *= (ms.zpbandflux(b) / normmagsys.zpbandflux(b))
    158 

/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/spectral.pyc in get_magsystem(name)
     37     if isinstance(name, MagSystem):
     38         return name
---> 39     return _MAGSYSTEMS.retrieve(name)
     40 
     41 

/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/_registry.pyc in retrieve(self, name, version)
    169             raise Exception(
    170                 "{0!r} not in registry. Registered names: '{1}'"
--> 171                 .format(name, "', '".join(registered_names)))
    172 
    173         # If version was specified but we don't have that specific version:

Exception: 'ab_b12' not in registry. Registered names: 'vega', 'ab', 'csp', 'bd17'

In [4]:
lc = sncosmo.read_lc('/Users/rbiswas/data/JLA_Calibration_OnlineData/cosmology/jla_light_curves/lc-SDSS14331.list',
                    format='salt2')

In [21]:
meta = lc.meta

In [23]:
lcp = lc.to_pandas()
#lcp['Filter'] = lcp['Filter'].apply(lambda x: 'sdss' + x[-1])
lcp['MagSys'] = 'ab'

In [24]:
lco  = LightCurve(lcp)

In [25]:
lcp.Filter.unique()


Out[25]:
array(['SDSS::g', 'SDSS::r', 'SDSS::i', 'SDSS::z', 'SDSS::u'], dtype=object)

In [26]:
bandDict = dict((x.lower(), x[:-3].lower() + x[-1]) for x in lcp.Filter.unique())

In [30]:
sncosmo.get_bandpass('mega_g')


Out[30]:
<Bandpass 'mega_g' at 0x1087a3fd0>

In [11]:
bandDict


Out[11]:
{'sdss::g': 'sdssg',
 'sdss::i': 'sdssi',
 'sdss::r': 'sdssr',
 'sdss::u': 'sdssu',
 'sdss::z': 'sdssz'}

In [12]:
LightCurve.remap_filters('SDSS::g', bandDict, True)


Out[12]:
'sdssg'

In [6]:
lco = LightCurve(lcp, bandNameDict=bandDict, ignore_case=True)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-f6912cc5add8> in <module>()
----> 1 lco = LightCurve(lcp, bandNameDict=bandDict, ignore_case=True)

NameError: name 'lcp' is not defined

In [14]:
lco.lightCurve.head()


Out[14]:
mjd flux fluxerr zp band zpsys
0 53974.332 -3.67 12.337 27.5 sdssg ab
1 53974.332 -3.77 18.608 27.5 sdssr ab
2 53974.332 -13.39 29.725 27.5 sdssi ab
3 53974.332 5.99 91.963 27.5 sdssz ab
4 53974.332 6.97 37.978 27.5 sdssu ab

In [15]:
dust = sncosmo.CCM89Dust()
model = sncosmo.Model(source='salt2', effects=[dust], effect_names=['mw'], effect_frames=['obs'])
model.set(mwebv=meta['MWEBV'])
#model.set(z=0.21)

In [16]:
import statsmodels

In [12]:
from scipy import stats

In [13]:
fitres = sncosmo.mcmc_lc(lco.snCosmoLC(), model=model, bounds=dict(z=(0.01,1.05), c=(-0.5, 0.5), x1=(-5.,5.)),
                         vparam_names=['x0', 'x1', 'c', 't0', 'z'],
                         priors=dict(z=lambda x: stats.norm.pdf(x, 1.13, 0.1)),
                         modelcov=True)


emcee: Exception while calling your likelihood function:
  params: [  5.30019214e-01   5.40109228e+04   4.62803751e-05   1.26256059e-04
  -1.91472324e-04]
  args: []
  kwargs: {}
  exception:
/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/fitting.py:132: RuntimeWarning: Dropping following bands from data: 'sdssz', 'sdssu', 'sdssg'(out of model wavelength range)
  "(out of model wavelength range)", RuntimeWarning)
Traceback (most recent call last):
  File "/usr/local/miniconda/lib/python2.7/site-packages/emcee/ensemble.py", line 519, in __call__
    return self.f(x, *self.args, **self.kwargs)
  File "/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/fitting.py", line 924, in lnprob
    logp += math.log(func(parameters[i]))
NameError: global name 'math' is not defined
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-ff4d1cfedd07> in <module>()
      2                          vparam_names=['x0', 'x1', 'c', 't0', 'z'],
      3                          priors=dict(z=lambda x: stats.norm.pdf(x, 1.13, 0.1)),
----> 4                          modelcov=True)

/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/fitting.pyc in mcmc_lc(data, model, vparam_names, bounds, priors, guess_amplitude, guess_t0, guess_z, minsnr, modelcov, nwalkers, nburn, nsamples, thin, a)
    942     # Run the sampler.
    943     sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, a=a)
--> 944     pos, prob, state = sampler.run_mcmc(pos, nburn)  # burn-in
    945     sampler.reset()
    946     sampler.run_mcmc(pos, nsamples, thin=thin)  # production run

/usr/local/miniconda/lib/python2.7/site-packages/emcee/sampler.pyc in run_mcmc(self, pos0, N, rstate0, lnprob0, **kwargs)
    170 
    171         for results in self.sample(pos0, lnprob0, rstate0, iterations=N,
--> 172                                    **kwargs):
    173             pass
    174 

/usr/local/miniconda/lib/python2.7/site-packages/emcee/ensemble.pyc in sample(self, p0, lnprob0, rstate0, blobs0, iterations, thin, storechain, mh_proposal)
    196         blobs = blobs0
    197         if lnprob is None:
--> 198             lnprob, blobs = self._get_lnprob(p)
    199 
    200         # Check to make sure that the probability function didn't return

/usr/local/miniconda/lib/python2.7/site-packages/emcee/ensemble.pyc in _get_lnprob(self, pos)
    380 
    381         # Run the log-probability calculations (optionally in parallel).
--> 382         results = list(M(self.lnprobfn, [p[i] for i in range(len(p))]))
    383 
    384         try:

/usr/local/miniconda/lib/python2.7/site-packages/emcee/ensemble.pyc in __call__(self, x)
    517     def __call__(self, x):
    518         try:
--> 519             return self.f(x, *self.args, **self.kwargs)
    520         except:
    521             import traceback

/Users/rbiswas/src/LSST/lsst_stack/DarwinX86/sncosmo/12.0-1-g5e70e90+2/lib/python/sncosmo-1.3.0-py2.7-macosx-10.6-x86_64.egg/sncosmo/fitting.pyc in lnprob(parameters)
    922 
    923         for i, func in idxpriors:
--> 924             logp += math.log(func(parameters[i]))
    925 
    926         return logp

NameError: global name 'math' is not defined

In [23]:
fitres = sncosmo.mcmc_lc(lco.snCosmoLC(), model=model, bounds=dict(z=(0.05, 0.8), c=(-0.5, 0.5), x1=(-5., 5.)),
                        vparam_names=['x0', 'x1', 'c', 't0', 'z'],
                        # priors=dict(z=lambda x: stats.norm.pdf(x,0.213, 0.2)), 
                        modelcov=True)

In [16]:
from analyzeSN import ResChar

In [17]:
mcmc_reschar = ResChar.fromSNCosmoRes(fitres)

In [18]:
mcmc_reschar.salt_samples().z.std()


Out[18]:
0.15658965129429681

In [19]:
np.histogram(mcmc_reschar.salt_samples().z)[1][np.histogram(mcmc_reschar.salt_samples().z)[0].argmax()]


Out[19]:
0.60148444622403763

In [20]:
fig, ax = plt.subplots()
mcmc_reschar.salt_samples().mu.hist(histtype='step',
                                    bins=np.arange(10., 14.,0.05),
                                    ax=ax)
#x = np.arange(0.1, 0.4, 0.001)
#ax.plot(x, 120. *stats.norm.pdf(x, 0.213, 0.4), 'r--')


Out[20]:
<matplotlib.axes._subplots.AxesSubplot at 0x10bc648d0>

In [21]:
fig, ax = plt.subplots()
mcmc_reschar.salt_samples().z.hist(histtype='step', 
                                   #bins=np.arange(0.126, 0.3126,0.005),
                                   ax=ax)
x = np.arange(0.1, 1.4, 0.001)
ax.plot(x, 120.*stats.norm.pdf(x, 1.13, 0.1), 'r--')


Out[21]:
[<matplotlib.lines.Line2D at 0x10bf53990>]

In [148]:
ax.plot(z, )


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-148-e4eaa140d4be> in <module>()
----> 1 ax.plot(z, )

NameError: name 'z' is not defined

In [ ]:
import seaborn as sns
sns.set_style('whitegrid')

In [149]:
sns.pairplot(data=mcmc_reschar.salt_samples())


Out[149]:
<seaborn.axisgrid.PairGrid at 0x10f650bd0>

In [131]:
ResChar.fromSNCosmoRes(fitres).parameters


Out[131]:
z            0.851218
t0       54021.692575
x0           0.000027
x1           3.054032
c           -0.352032
mwebv        0.022100
mwr_v        3.100000
dtype: float64

In [132]:
fig = sncosmo.plot_lc(lco.snCosmoLC(), model=fitres[1])



In [15]:
fig


Out[15]:

Setting model parameters from JLA_paramslc.txt by hand


In [16]:
model.set(x1=0.141351)
model.set(z=0.213)
model.set(c=0.0142)

In [17]:
print model


<Model at 0x10a48b410>
source:
  class      : SALT2Source
  name       : 'salt2'
  version    : 2.4
  phases     : [-20, .., 50] days
  wavelengths: [2000, .., 9200] Angstroms
effect (name='mw' frame='obs'):
  class           : CCM89Dust
  wavelength range: [1000, 33333.3] Angstroms
parameters:
  z     = 0.21299999999999999
  t0    = 0.0
  x0    = 1.0
  x1    = 0.141351
  c     = 0.014200000000000001
  mwebv = 0.022100000000000002
  mwr_v = 3.1000000000000001

In [20]:
fitres_JLAparams = sncosmo.fit_lc(lco.snCosmoLC(), model=model, vparam_names=['x0', 't0'], modelcov=True)

In [21]:
fitres_JLAparams[0].parameters


Out[21]:
array([  2.13000000e-01,   5.40134824e+04,   5.79360360e-05,
         1.41351000e-01,   1.42000000e-02,   2.21000000e-02,
         3.10000000e+00])

In [22]:
fig = sncosmo.plot_lc(lco.snCosmoLC(), model=(fitres[1], fitres_JLAparams[1]))

In [23]:
fig


Out[23]:

In [20]:
mcmcres_JLAparams = sncosmo.mcmc_lc(lco.snCosmoLC(), model=model, vparam_names=['x0', 't0'], modelcov=True)

In [21]:
ResChar.fromSNCosmoRes(mcmcres_JLAparams).parameters


Out[21]:
z            0.213000
t0       54013.470555
x0           0.000058
x1           0.141351
c            0.014200
mwebv        0.022100
mwr_v        3.100000
dtype: float64

In [22]:
sncosmo.plot_lc(lco.snCosmoLC(), model=[fitres[1], fitres_JLAparams[1], mcmcres_JLAparams[1]], color='k')


Out[22]:

In [27]:
fitres_JLAparams = sncosmo.fit_lc(lco.snCosmoLC(), model=model, vparam_names=['x0', 't0'], modelcov=False, minsnr=5.)

In [28]:
fitres_JLAparams[0].parameters


Out[28]:
array([  2.13000000e-01,   5.40134475e+04,   5.83488019e-05,
         1.41351000e-01,   1.42000000e-02,   2.21000000e-02,
         3.10000000e+00])

In [ ]: