In [13]:
from astropy.io import fits
import pandas as pd
import numpy as np
from sqlalchemy import create_engine
%matplotlib inline
import matplotlib.pyplot as plt
from astropy.table import Table

In [44]:
import sncosmo

In [2]:
import analyzeSN as ans

In [3]:
SDSS_summary_file = ans.snanaSims.SnanaSims.snanadatafile(snanafileroot='SDSS_allCandidates+BOSS', location='/Users/rbiswas/data/SNDATA/lcmerge')

In [19]:
summary = Table(fits.open(SDSS_summary_file)[1].data)

In [66]:
SNIDs = map (lambda x: x.strip(), summary[summary['SNTYPE'] == 120]['SNID'])

In [33]:
sdss = ans.snanaSims.SnanaSims.fromSNANAfileroot(snanafileroot='SDSS_allCandidates+BOSS', location='/Users/rbiswas/data/SNDATA/lcmerge'
                                                 )


/Users/rbiswas/data/SNDATA/lcmerge/SDSS_allCandidates+BOSS_HEAD.FITS

In [42]:
sne = ans.snanaSims.SnanaSims.fromSNANAfileroot(snanafileroot='SDSS_allCandidates+BOSS', 
                                                location='/Users/rbiswas/data/SNDATA/lcmerge',
                                                snids=SNIDs)


/Users/rbiswas/data/SNDATA/lcmerge/SDSS_allCandidates+BOSS_HEAD.FITS

In [43]:
sn0 = ans.snanaSims.SnanaSims.reformat_SNANASN(sne.snList[0], snanaBands=['u','g','r','i','z'], replacements=['sdssu','sdssg', 'sdssr', 'sdssi', 'sdssz'])

In [108]:
sn1 = ans.snanaSims.SnanaSims.reformat_SNANASN(sne.snList[1], snanaBands=['u','g','r','i','z'], replacements=['sdssu','sdssg', 'sdssr', 'sdssi', 'sdssz'])

In [72]:
sn0[:5]


Out[72]:
MJDFIELDTELESCOPEPHOTFLAGfluxfluxerrMAGMAGERRPSF_SIG1PSF_SIG2PSF_RATIOSKY_SIGRDNOISEZEROPTZEROPT_ERRGAINZPZPSYSband
53616.316406282NSDSS11562721.42233.39618.9130.0931.4332.9180.0934.274.25327.120.014.0527.5absdssg
53616.316406282NSDSS11643577.63442.67818.6160.1341.4133.250.075.964.25225.3550.0014.72527.5absdssr
53616.316406282NSDSS11562436.55254.07519.0330.1131.1272.6510.0857.212.99225.5630.0024.6427.5absdssi
53616.316406282NSDSS361143.38333.67419.850.311.3083.4080.0477.544.69825.4270.0063.4827.5absdssz
53616.316406282NSDSS1156611.75180.76120.5330.3161.4483.0650.0825.144.48326.0990.01.4727.5absdssu

In [70]:
fig_data = sncosmo.plot_lc(sn0);



In [62]:
fig_data.savefig('SN_data.pdf')

In [63]:
print sn0.meta


OrderedDict([('SNID', '722'), ('IAUC', '2005ed  '), ('FAKE', 0), ('RA', 0.70565999999999995), ('DECL', 0.75116300000000003), ('PIXSIZE', 0.40000001), ('SNTYPE', 120), ('NOBS', 105), ('PTROBS_MIN', 3824), ('PTROBS_MAX', 3928), ('MWEBV', 0.026000001), ('MWEBV_ERR', 0.0043000001), ('REDSHIFT_HELIO', 0.086341001), ('REDSHIFT_HELIO_ERR', 3.6000001e-05), ('REDSHIFT_FINAL', 0.085037999), ('REDSHIFT_FINAL_ERR', 3.6000001e-05), ('HOSTGAL_OBJID', 1237657191979286528), ('HOSTGAL_PHOTOZ', 0.087899998), ('HOSTGAL_PHOTOZ_ERR', -9.0), ('HOSTGAL_SNSEP', 2.2219999), ('HOSTGAL_LOGMASS', 0.0), ('HOSTGAL_LOGMASS_ERR', 0.0), ('PEAKMJD', 53609.0), ('SEARCH_TYPE', 120)])

In [141]:
modelz = sncosmo.Model(source='salt2')
modelz.set(z = sn0.meta['REDSHIFT_FINAL'])
fitted_resz, fitted_modelz = sncosmo.fit_lc(sn0, modelz, vparam_names=['t0', 'x0', 'x1', 'c'])
fig_modelz = sncosmo.plot_lc(sn0, fitted_modelz, show_pulls=True)



In [145]:
fitted_resz.parameters


Out[145]:
array([  8.50379989e-02,   5.36134133e+04,   4.63145275e-04,
        -9.16142353e-01,  -1.13496200e-02])

In [82]:
fig_modelz.savefig('SN_data_fit.pdf')

In [83]:
fitted_resz.errors


Out[83]:
OrderedDict([('t0', 0.27924119317685836), ('x0', 1.525427341184638e-05), ('x1', 0.11545855720456534), ('c', 0.017422062840144397)])

In [74]:
mcmc_resz


Out[74]:
             vparam_names: ['t0', 'x0', 'x1', 'c']
                   errors: OrderedDict([('t0', 0.35116271215957978), ('x0', 1.8135242844808131e-05), ('x1', 0.11898216947736397), ('c', 0.01967591634566334)])
                  samples: array([[  5.36145164e+04,   4.10110601e-04,  -4.96340289e-01,
          2.96097592e-02],
       [  5.36145164e+04,   4.10110601e-04,  -4.96340289e-01,
          2.96097592e-02],
       [  5.36145209e+04,   4.10368620e-04,  -5.43051894e-01,
          2.99999252e-02],
       ..., 
       [  5.36132378e+04,   4.69158332e-04,  -9.38790605e-01,
         -1.54022249e-02],
       [  5.36130909e+04,   4.66246996e-04,  -9.31236998e-01,
         -1.21236507e-02],
       [  5.36130909e+04,   4.66246996e-04,  -9.31236998e-01,
         -1.21236507e-02]])
               parameters: array([  8.50379989e-02,   5.36134801e+04,   4.60586372e-04,
        -9.03485239e-01,  -9.53302887e-03])
              param_names: ['z', 't0', 'x0', 'x1', 'c']
 mean_acceptance_fraction: 0.59369999999999989
               covariance: array([[  1.23315250e-01,  -5.36054029e-06,   8.81914128e-03,
          4.93441153e-03],
       [ -5.36054029e-06,   3.28887033e-10,  -1.31231813e-06,
         -3.29665651e-07],
       [  8.81914128e-03,  -1.31231813e-06,   1.41567567e-02,
          1.24888810e-03],
       [  4.93441153e-03,  -3.29665651e-07,   1.24888810e-03,
          3.87141684e-04]])

In [73]:
mcmc_resz, mcmc_modelz = sncosmo.mcmc_lc(sn0, modelz, vparam_names=['t0', 'x0', 'x1', 'c'])

In [75]:
import triangle

In [159]:
figure_mcmc = triangle.corner(mcmc_resz.samples, labels=mcmc_resz.vparam_names, range=4*[0.95],
                         show_titles=True, title_args={"fontsize": 12})

figure_mcmc.gca().annotate("mcmc sampling", xy=(0.5, 1.0), xycoords="figure fraction",
                      xytext=(0, -5), textcoords="offset points",
                      ha="center", va="top")

axes = figure_mcmc.axes



In [84]:
figure_mcmc.savefig('mcmc_fig.pdf')

In [80]:
fig_comp = sncosmo.plot_lc(sn0, model=[fitted_modelz, mcmc_modelz], model_label=['MaxL', 'MCMC'])



In [85]:
fig_comp.savefig('fig_comp.pdf')

In [147]:
nest_res, nest_model = sncosmo.nest_lc(sn0, modelz, vparam_names=['t0', 'x0', 'x1', 'c'], 
                                    bounds={'c':(-0.3, 0.3), 'x1':(-3.0, 3.0)}, guess_amplitude_bound=True, minsnr=3.0, 
                                    verbose=True)


 iter=  2431 logz=-138.156004calls=6558 time= 14.312s

In [162]:
figure_nest = triangle.corner(nest_res.samples, weights=nest_res.weights, labels=nest_res.vparam_names, range=4*[0.99],
                              bins=20)

# figure_nest.gca().annotate("nest sampling", xy=(0.5, 1.0), xycoords="figure fraction",
#                      xytext=(0, -5), textcoords="offset points",
#                      ha="center", va="top")



In [156]:
triangle.corner?

In [97]:
model = sncosmo.Model(source='salt2')

In [99]:
fitted_res, fitted_model = sncosmo.fit_lc(sn0, model, vparam_names=['z', 't0', 'x0', 'x1', 'c'], bounds={'z':(0., 1.2)})
fig_model = sncosmo.plot_lc(sn0, fitted_model, show_pulls=True)


/Users/rbiswas/.local/lib/python2.7/site-packages/sncosmo-1.1.dev516-py2.7-macosx-10.5-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)

In [100]:
print fitted_model


<Model at 0x10b6714d0>
source:
  class      : SALT2Source
  name       : 'salt2'
  version    : 2.4
  phases     : [-20, .., 50] days
  wavelengths: [2000, .., 9200] Angstroms
parameters:
  z  = 0.50537863412860085
  t0 = 53613.082380537286
  x0 = 0.00020210570681280184
  x1 = 6.1862047718297601
  c  = -0.3850243045219397

In [101]:
mcmc_res, mcmc_model = sncosmo.mcmc_lc(sn0, model, vparam_names=['z', 't0', 'x0', 'x1', 'c'], bounds={'z':(0., 1.2)})
fig_mcmc_model = sncosmo.plot_lc(sn0, mcmc_model, show_pulls=True)



In [102]:
print mcmc_model


<Model at 0x10d54d850>
source:
  class      : SALT2Source
  name       : 'salt2'
  version    : 2.4
  phases     : [-20, .., 50] days
  wavelengths: [2000, .., 9200] Angstroms
parameters:
  z  = 0.63218577324082303
  t0 = 53616.273191568042
  x0 = 0.0001574520167012656
  x1 = 6.3340704677295241
  c  = -0.4413991859321762

In [104]:
mcmc_res


Out[104]:
             vparam_names: ['z', 't0', 'x0', 'x1', 'c']
                   errors: OrderedDict([('z', 0.021270960051247321), ('t0', 0.72396832720842397), ('x0', 8.0933250635547245e-06), ('x1', 0.69552657232405724), ('c', 0.060384558293747108)])
                  samples: array([[  5.94833245e-01,   5.36166066e+04,   1.65097651e-04,
          5.94610338e+00,  -3.63502853e-01],
       [  5.94733913e-01,   5.36165721e+04,   1.65754881e-04,
          5.91543318e+00,  -3.70029981e-01],
       [  5.94728620e-01,   5.36165695e+04,   1.65888709e-04,
          5.90656242e+00,  -3.71065119e-01],
       ..., 
       [  6.49638656e-01,   5.36160125e+04,   1.51445952e-04,
          6.70700546e+00,  -4.85178673e-01],
       [  6.49882706e-01,   5.36160490e+04,   1.52519486e-04,
          6.55339997e+00,  -4.99952852e-01],
       [  6.49882706e-01,   5.36160490e+04,   1.52519486e-04,
          6.55339997e+00,  -4.99952852e-01]])
               parameters: array([  6.32185773e-01,   5.36162732e+04,   1.57452017e-04,
         6.33407047e+00,  -4.41399186e-01])
              param_names: ['z', 't0', 'x0', 'x1', 'c']
 mean_acceptance_fraction: 0.55300000000000016
               covariance: array([[  4.52453742e-04,   6.13347168e-03,  -1.35183356e-07,
          2.92940011e-03,  -1.10630070e-03],
       [  6.13347168e-03,   5.24130139e-01,  -6.53769420e-07,
         -2.84005519e-01,  -2.23569401e-02],
       [ -1.35183356e-07,  -6.53769420e-07,   6.55019106e-11,
         -3.75188073e-06,   2.90709361e-07],
       [  2.92940011e-03,  -2.84005519e-01,  -3.75188073e-06,
          4.83757213e-01,   4.05177895e-03],
       [ -1.10630070e-03,  -2.23569401e-02,   2.90709361e-07,
          4.05177895e-03,   3.64629488e-03]])

In [105]:
figure_mcmc = triangle.corner(mcmc_res.samples, labels=mcmc_res.vparam_names,
                         show_titles=True, title_args={"fontsize": 12})

figure_mcmc.gca().annotate("mcmc sampling", xy=(0.5, 1.0), xycoords="figure fraction",
                      xytext=(0, -5), textcoords="offset points",
                      ha="center", va="top")

axes = figure_mcmc.axes



In [106]:
nest_res, nest_model = sncosmo.mcmc_lc(sn0, model, vparam_names=['z', 't0', 'x0', 'x1', 'c'], bounds={'z':(0., 1.2)})
fig_nest_model = sncosmo.plot_lc(sn0, nest_model, show_pulls=True)



In [113]:
nest_res_small, nest_model_small = sncosmo.mcmc_lc(sn0, model, vparam_names=['z', 't0', 'x0', 'x1', 'c'], bounds={'z':(0., 0.3)})
fig_nest_model_small = sncosmo.plot_lc(sn0, nest_model_small, show_pulls=True)



In [140]:
fitted_resz.parameters


Out[140]:
array([  8.50379989e-02,   5.36134133e+04,   4.63145275e-04,
        -9.16142353e-01,  -1.13496200e-02])

In [138]:
nest_res_small.parameters


Out[138]:
array([  1.13283425e-01,   5.36127746e+04,   4.67917577e-04,
        -6.25916911e-01,  -5.47341339e-02])

In [136]:
nest_res_small.errors


Out[136]:
OrderedDict([('z', 0.032940945706096181), ('t0', 0.88114021982012847), ('x0', 2.3371803237338381e-05), ('x1', 0.3649718135313233), ('c', 0.061031847496703068)])

In [130]:
#sn1.rename_column('ZP', 'zp')
sn1.rename_column('ZPSYS', 'zpsys')

In [131]:
sn1


Out[131]:
MJDFIELDTELESCOPEPHOTFLAGfluxfluxerrMAGMAGERRPSF_SIG1PSF_SIG2PSF_RATIOSKY_SIGRDNOISEZEROPTZEROPT_ERRGAINzpzpsysband
53616.355468882NSDSS01505.1736.0519.5560.0261.6683.4190.1454.084.25328.3730.0034.0527.5absdssg
53616.355468882NSDSS01732.9238.31319.4030.0241.5323.5050.0885.294.25228.1880.0034.72527.5absdssr
53616.355468882NSDSS01228.8749.83619.7760.0441.4673.6790.0656.8812.99227.8560.0034.6427.5absdssi
53616.355468882NSDSS32562.16157.12320.6070.291.5063.5660.0697.114.69826.5330.0053.4827.5absdssz
53616.355468882NSDSS8178.3849.32221.8650.2931.6853.4730.1454.924.48327.4470.0261.4727.5absdssu
53626.332031282NSDSS0532.4520.12320.6840.0411.2962.9760.0594.374.25328.4370.0034.0527.5absdssg
53626.332031282NSDSS01014.727.11419.9840.0291.0552.4470.0815.294.25228.2330.0054.72527.5absdssr
53626.332031282NSDSS0820.7140.11720.2140.0531.0212.5330.0696.6612.99227.8870.0024.6427.5absdssi
53626.332031282NSDSS32474.46128.1820.7840.2771.3373.6190.0426.434.69826.5230.013.4827.5absdssz
53626.332031282NSDSS3255.1644.65423.0820.7281.3673.0790.0784.84.48327.5840.0161.4727.5absdssu
53628.316406282NSDSS1369.2743.68721.0810.1281.4633.5440.05112.94.25328.3720.0074.0527.5absdssg
.........................................................
53698.269531282NSDSS8-11.6150.04225.0731.3851.4683.310.0644.584.48327.40180.02351.4727.5absdssu
53699.292968882NSDSS3261.9317.92522.9980.2981.1012.2440.1674.084.25328.3790.0054.0527.5absdssg
53699.292968882NSDSS3275.2625.44922.7820.3441.1872.7320.0515.554.25228.2020.0034.72527.5absdssr
53699.292968882NSDSS3252.7941.20923.0850.6591.0652.5180.0637.2512.99227.9010.0044.6427.5absdssi
53699.292968882NSDSS0-75.62109.56223.360.6720.9662.4510.0627.44.69826.6340.0063.4827.5absdssz
53699.292968882NSDSS085.8448.00222.6380.5531.2892.6230.0934.424.48327.4050.0071.4727.5absdssu
53704.257812582NSDSS3289.322.02422.6120.261.5573.1260.0914.054.25328.4020.0024.0527.5absdssg
53704.257812582NSDSS3290.8328.81422.5860.3281.4993.3870.065.434.25228.2190.0014.72527.5absdssr
53704.257812582NSDSS3267.1349.62522.8620.6651.4013.5930.0457.9612.99227.9010.0024.6427.5absdssi
53704.257812582NSDSS40146.46143.03921.8780.6971.3443.410.0518.074.69826.6230.0033.4827.5absdssz
53704.257812582NSDSS8-45.9455.21726.0140.9771.6663.4110.0814.54.48327.5310.0021.4727.5absdssu

In [111]:
nest_res, nest_model = sncosmo.mcmc_lc(sn1, model, vparam_names=['z', 't0', 'x0', 'x1', 'c'], bounds={'z':(0., 0.3)})
fig_nest_model = sncosmo.plot_lc(sn1, nest_model, show_pulls=True)



In [125]:
sn1['zp']


---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-125-165112bb813f> in <module>()
----> 1 sn1['zp']

/usr/local/manual/anaconda/lib/python2.7/site-packages/astropy/table/table.pyc in __getitem__(self, item)
    746     def __getitem__(self, item):
    747         if isinstance(item, six.string_types):
--> 748             return self.columns[item]
    749         elif isinstance(item, (int, np.integer)):
    750             return self.Row(self, item)

/usr/local/manual/anaconda/lib/python2.7/site-packages/astropy/table/table.pyc in __getitem__(self, item)
     77         """
     78         if isinstance(item, six.string_types):
---> 79             return OrderedDict.__getitem__(self, item)
     80         elif isinstance(item, int):
     81             return self.values()[item]

KeyError: 'zp'

In [132]:
sncosmo.chisq(sn1, nest_model)


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-132-6c79d718f203> in <module>()
----> 1 sncosmo.chisq(sn1, nest_model)

/Users/rbiswas/.local/lib/python2.7/site-packages/sncosmo-1.1.dev516-py2.7-macosx-10.5-x86_64.egg/sncosmo/fitting.pyc in chisq(data, model, modelcov)
     63     """
     64     data = standardize_data(data)
---> 65     return _chisq(data, model, modelcov=modelcov)
     66 
     67 

/Users/rbiswas/.local/lib/python2.7/site-packages/sncosmo-1.1.dev516-py2.7-macosx-10.5-x86_64.egg/sncosmo/fitting.pyc in _chisq(data, model, modelcov)
     39     else:
     40         mflux = model.bandflux(data['band'], data['time'],
---> 41                                zp=data['zp'], zpsys=data['zpsys'])
     42         return np.sum(((data['flux'] - mflux) / data['fluxerr'])**2)
     43 

/Users/rbiswas/.local/lib/python2.7/site-packages/sncosmo-1.1.dev516-py2.7-macosx-10.5-x86_64.egg/sncosmo/models.pyc in bandflux(self, band, time, zp, zpsys)
   1293         except ValueError as e:
   1294             _check_for_fitpack_error(e, time, 'time')
-> 1295             raise e
   1296 
   1297     def _bandflux_rcov(self, band, time):

ValueError: bandpass 'sdssu' [2940, .., 7940] outside spectral range [2953.85, .., 13587.7]

In [118]:
type(sn1)


Out[118]:
astropy.table.table.Table

In [119]:
from astropy.table import Table

In [121]:
isinstance(sn1, Table)


Out[121]:
True

In [123]:
help(sncosmo.chisq)


Help on function chisq in module sncosmo.fitting:

chisq(data, model, modelcov=False)
    Calculate chisq statistic for the model, given the data.
    
    Parameters
    ----------
    model : `~sncosmo.Model`
    data : `~astropy.table.Table` or `~numpy.ndarray` or `dict`
        Table of photometric data. Must include certain columns.
        See the "Photometric Data" section of the documentation for
        required columns.
    modelcov : bool
        Include model covariance? Calls ``model.bandfluxcov`` method
        instead of ``model.bandflux``. The source in the model must therefore
        implement covariance.
    
    Returns
    -------
    chisq : float


In [151]:
len(nest_res.vparam_names)


Out[151]:
4

In [152]:
nest_res


Out[152]:
        niter: 2432
       errors: OrderedDict([('t0', 0.2872317471576008), ('x0', 1.5693132174788232e-05), ('x1', 0.11346429564081159), ('c', 0.017725694665084411)])
   param_dict: OrderedDict([('z', 0.085037998855113983), ('t0', 53613.409261921457), ('x0', 0.00046345059029101102), ('x1', -0.91566806594514438), ('c', -0.011487138383286079)])
   parameters: array([  8.50379989e-02,   5.36134093e+04,   4.63450590e-04,
        -9.15668066e-01,  -1.14871384e-02])
       bounds: {'c': (-0.3, 0.3), 'x1': (-3.0, 3.0), 'x0': (0.0, 0.0057901827652364762), 't0': (53562.064506307244, 53725.919509977102)}
            h: 16.461970749754641
     logprior: array([ -4.61016602,  -4.62016602,  -4.63016602, ..., -29.60517019,
       -29.60517019, -29.60517019])
         logz: -138.1545109966859
         ndof: 80
        ncall: 6558
 vparam_names: ['t0', 'x0', 'x1', 'c']
      weights: array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
         1.42910351e-05,   1.53752326e-05,   1.45962895e-05])
      samples: array([[  5.36466805e+04,   5.69231560e-03,   2.55536102e+00,
         -1.33199768e-01],
       [  5.36780593e+04,   4.09178797e-03,   1.48611524e+00,
         -8.22127902e-02],
       [  5.36750274e+04,   3.54027363e-03,   1.48462324e+00,
         -1.24438691e-01],
       ..., 
       [  5.36134752e+04,   4.58570734e-04,  -9.09544729e-01,
         -7.13654288e-03],
       [  5.36134640e+04,   4.61024519e-04,  -9.17165931e-01,
         -9.05829250e-03],
       [  5.36135077e+04,   4.57965369e-04,  -9.03513110e-01,
         -6.43367263e-03]])
         time: 14.312282800674438
      logzerr: 0.4057335424851468
   covariance: array([[  8.25020766e-02,  -3.40439025e-06,   4.98050528e-03,
          3.07205771e-03],
       [ -3.40439025e-06,   2.46274397e-10,  -1.11621020e-06,
         -2.53338202e-07],
       [  4.98050528e-03,  -1.11621020e-06,   1.28741464e-02,
          1.07616558e-03],
       [  3.07205771e-03,  -2.53338202e-07,   1.07616558e-03,
          3.14200251e-04]])
         logl: array([ -9.63308642e+06,  -8.37854954e+06,  -6.52095339e+06, ...,
        -1.19705219e+02,  -1.19632093e+02,  -1.19684084e+02])

In [ ]: