In [1]:

    

%matplotlib inline

import collections
import pandas as pd
import sncosmo
import triangle

import lc


    

In [2]:

    

import matplotlib.pyplot as plt
import os
import numpy as np
from astropy.units import Unit
from astropy.table import Table


    

In [3]:

    

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 [23]:

    

sncosmoband


    

    
Out[23]:





<Bandpass 'LSST_y' at 0x10546c350>

In [6]:

    

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)

data = sncosmo.read_lc('sn.dat')
vparams = ['t0', 'x0', 'x1', 'c']


    

In [9]:

    

data.meta


    

    
Out[9]:





OrderedDict([('SNID', 44), ('IAUC', 'NULL'), ('FAKE', 2), ('RA', 209.882553), ('DECL', -12.400555), ('PIXSIZE', 0.2), ('NXPIX', -9), ('NYPIX', -9), ('SNTYPE', 1), ('NOBS', 83), ('PTROBS_MIN', 73), ('PTROBS_MAX', 155), ('MWEBV', 0.0829484), ('MWEBV_ERR', 0.0132717), ('REDSHIFT_HELIO', 1.06465), ('REDSHIFT_HELIO_ERR', 0.0005), ('REDSHIFT_FINAL', 1.06656), ('REDSHIFT_FINAL_ERR', 0.0005), ('HOSTGAL_OBJID', 29144), ('HOSTGAL_PHOTOZ', 1.20071), ('HOSTGAL_PHOTOZ_ERR', 0.05), ('HOSTGAL_SPECZ', 0.0), ('HOSTGAL_SPECZ_ERR', 0.0), ('HOSTGAL_SNSEP', -999.0), ('HOSTGAL_LOGMASS', -9.0), ('HOSTGAL_LOGMASS_ERR', -9.0), ('HOSTGAL_MAG_u', 999.0), ('HOSTGAL_MAG_g', 999.0), ('HOSTGAL_MAG_r', 999.0), ('HOSTGAL_MAG_i', 999.0), ('HOSTGAL_MAG_z', 999.0), ('HOSTGAL_MAG_Y', 999.0), ('HOSTGAL_SB_FLUXCAL_u', -9.0), ('HOSTGAL_SB_FLUXCAL_g', -9.0), ('HOSTGAL_SB_FLUXCAL_r', -9.0), ('HOSTGAL_SB_FLUXCAL_i', -9.0), ('HOSTGAL_SB_FLUXCAL_z', -9.0), ('HOSTGAL_SB_FLUXCAL_Y', -9.0), ('PEAKMJD', 50915.8), ('SEARCH_TYPE', -9), ('SIM_MODEL_NAME', 'SALT2.Guy10_LAMOPEN'), ('SIM_MODEL_INDEX', 6), ('SIM_TYPE_INDEX', 1), ('SIM_TYPE_NAME', 'Ia'), ('SIM_NON1a', 0), ('SIM_LIBID', 2082), ('SIM_SEARCHEFF_MASK', 3), ('SIM_REDSHIFT_HELIO', 1.06371), ('SIM_REDSHIFT_CMB', 1.06562), ('SIM_VPEC', -0.0), ('SIM_DLMU', 44.2709), ('SIM_RA', 209.882553101), ('SIM_DECL', -12.400554657), ('SIM_MWEBV', 0.0821098), ('SIM_PEAKMJD', 50916.5), ('SIM_SALT2x0', 1.695e-06), ('SIM_SALT2x1', -0.708466), ('SIM_SALT2c', 0.0178018), ('SIM_SALT2mB', 25.0621), ('SIM_SALT2alpha', 0.14), ('SIM_SALT2beta', 3.2), ('SIM_PEAKMAG_u', -9.0), ('SIM_PEAKMAG_g', -9.0), ('SIM_PEAKMAG_r', 26.5509), ('SIM_PEAKMAG_i', 25.0618), ('SIM_PEAKMAG_z', 24.534), ('SIM_PEAKMAG_Y', 24.5804), ('SIM_EXPOSURE_u', 1.0), ('SIM_EXPOSURE_g', 1.0), ('SIM_EXPOSURE_r', 1.0), ('SIM_EXPOSURE_i', 1.0), ('SIM_EXPOSURE_z', 1.0), ('SIM_EXPOSURE_Y', 1.0), ('SIM_GALFRAC_u', -9.0), ('SIM_GALFRAC_g', -9.0), ('SIM_GALFRAC_r', -9.0), ('SIM_GALFRAC_i', -9.0), ('SIM_GALFRAC_z', -9.0), ('SIM_GALFRAC_Y', -9.0)])

In [10]:

    

test = lc.LC(model, data, vparams)
iotest = lc.LC(model, data, vparams)


    

In [11]:

    

test.metadata(), iotest.metadata()


    

    




<Model at 0x1099f2b10>
source:
  class      : SALT2Source
  name       : 'salt2-extended'
  version    : 1.0
  phases     : [-20, .., 50] days
  wavelengths: [300, .., 18000] Angstroms
parameters:
  z  = 1.0637099999999999
  t0 = 50916.5
  x0 = 1.6950000000000001e-06
  x1 = -0.70846600000000004
  c  = 0.0178018
<Model at 0x1099f2b10>
source:
  class      : SALT2Source
  name       : 'salt2-extended'
  version    : 1.0
  phases     : [-20, .., 50] days
  wavelengths: [300, .., 18000] Angstroms
parameters:
  z  = 1.0637099999999999
  t0 = 50916.5
  x0 = 1.6950000000000001e-06
  x1 = -0.70846600000000004
  c  = 0.0178018






    
Out[11]:





(None, None)

In [12]:

    

test.fitOut


    

    




running chi^2 fit






    
Out[12]:





(       errors: OrderedDict([('t0', 0.869705464632716), ('x0', 6.352280833289825e-08), ('x1', 0.5247654629885447), ('c', 0.048283506458290146)])
   parameters: array([  1.06371000e+00,   5.09155364e+04,   1.19022208e-06,
        -1.37443229e+00,  -2.74520718e-02])
      success: True
         ndof: 79
   covariance: array([[  7.56396247e-01,  -1.47415535e-08,   9.27513342e-02,
          8.97227286e-03],
       [ -1.47415535e-08,   4.03514718e-15,  -1.57647791e-08,
          8.09694619e-10],
       [  9.27513342e-02,  -1.57647791e-08,   2.79008873e-01,
          1.00106837e-02],
       [  8.97227286e-03,   8.09694619e-10,   1.00106837e-02,
          2.35188405e-03]])
 vparam_names: ['t0', 'x0', 'x1', 'c']
        chisq: 94.67428083189289
  param_names: ['z', 't0', 'x0', 'x1', 'c']
      message: 'Minimization exited successfully.'
        ncall: 98,
 <sncosmo.models.Model at 0x103811610>)

In [8]:

    

test.mcmcOut


    

    




running MCMC fit






    
Out[8]:





(             vparam_names: ['t0', 'x0', 'x1', 'c']
                   errors: OrderedDict([('t0', 0.84052291965928672), ('x0', 6.6097763429030765e-08), ('x1', 0.54376647284091717), ('c', 0.051716035384276604)])
                  samples: array([[  5.09172430e+04,   1.15286766e-06,  -5.73551900e-01,
          7.80925076e-02],
       [  5.09169532e+04,   1.15937672e-06,  -6.53440242e-01,
          6.60318272e-02],
       [  5.09169532e+04,   1.15937672e-06,  -6.53440242e-01,
          6.60318272e-02],
       ..., 
       [  5.09158980e+04,   1.22719862e-06,  -1.87470109e+00,
         -1.68128260e-02],
       [  5.09155473e+04,   1.24432644e-06,  -1.67843172e+00,
          4.68350001e-03],
       [  5.09155473e+04,   1.24432644e-06,  -1.67843172e+00,
          4.68350001e-03]])
               parameters: array([  1.06371000e+00,   5.09154497e+04,   1.19142568e-06,
        -1.26646862e+00,  -1.10709753e-02])
              param_names: ['z', 't0', 'x0', 'x1', 'c']
 mean_acceptance_fraction: 0.59349999999999992
               covariance: array([[  7.06478778e-01,  -1.09831320e-08,   1.11635706e-01,
          9.94109082e-03],
       [ -1.09831320e-08,   4.36891433e-15,  -1.94720479e-08,
          9.84437375e-10],
       [  1.11635706e-01,  -1.94720479e-08,   2.95681977e-01,
          1.01745451e-02],
       [  9.94109082e-03,   9.84437375e-10,   1.01745451e-02,
          2.67454832e-03]]),
 <sncosmo.models.Model at 0x109ac65d0>)

In [9]:

    

test.nestOut


    

    




running nest fit
 iter=  1805 logz=-60.447548calls=6163 time=398.265s






    
Out[9]:





(        niter: 1806
       errors: OrderedDict([('t0', 0.86438314449960862), ('x0', 6.197342094873836e-08), ('x1', 0.54924416000489951), ('c', 0.04960572640044595)])
   param_dict: OrderedDict([('z', 1.0637099999999999), ('t0', 50915.602091922599), ('x0', 1.1863274802416421e-06), ('x1', -1.2611093831758817), ('c', -0.013469999689044585)])
   parameters: array([  1.06371000e+00,   5.09156021e+04,   1.18632748e-06,
        -1.26110938e+00,  -1.34699997e-02])
       bounds: {'c': (-0.3, 0.3), 'x1': (-3.0, 3.0), 'x0': (0.0, 1.4249497311936984e-05), 't0': (50774.025500000011, 51071.31319999999)}
            h: 11.062512885719315
     logprior: array([ -4.61016602,  -4.62016602,  -4.63016602, ..., -22.66517019,
       -22.66517019, -22.66517019])
         logz: -60.441002229837324
         ndof: 79
        ncall: 6163
 vparam_names: ['t0', 'x0', 'x1', 'c']
      weights: array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
         6.27371338e-05,   6.52131300e-05,   6.57808167e-05])
      samples: array([[  5.09150267e+04,   1.36766909e-05,   8.71338960e-01,
          2.42580968e-02],
       [  5.09108833e+04,   1.31580838e-05,  -1.04616311e+00,
         -1.34483275e-01],
       [  5.09755581e+04,   1.21386874e-05,   2.39309984e+00,
         -1.57536244e-01],
       ..., 
       [  5.09159287e+04,   1.18173883e-06,  -1.23987817e+00,
         -1.59102013e-02],
       [  5.09153033e+04,   1.17754694e-06,  -1.31047116e+00,
         -3.32303638e-02],
       [  5.09155558e+04,   1.19905724e-06,  -1.55125633e+00,
         -3.69537136e-02]])
         time: 398.26488494873047
      logzerr: 0.33260356110119016
   covariance: array([[  7.47158220e-01,  -1.22990009e-08,   6.04041064e-02,
          9.30440893e-03],
       [ -1.22990009e-08,   3.84070490e-15,  -1.49347510e-08,
          7.36404831e-10],
       [  6.04041064e-02,  -1.49347510e-08,   3.01669147e-01,
          1.04489323e-02],
       [  9.30440893e-03,   7.36404831e-10,   1.04489323e-02,
          2.46072809e-03]])
         logl: array([ -5.04238964e+04,  -3.96944503e+04,  -3.50742919e+04, ...,
        -4.74523891e+01,  -4.74136818e+01,  -4.74050143e+01]),
 <sncosmo.models.Model at 0x103921dd0>)

In [11]:

    

test.writeFits('fits.hdf5')


    

In [12]:

    

iotest.fitOut, iotest.mcmcOut, iotest.nestOut = iotest.readFits('fits.hdf5')


    

In [13]:

    

iotest.fitOut


    

    
Out[13]:





(       errors: OrderedDict([('t0', 0.86970546463271603), ('x0', 6.3522808332898253e-08), ('x1', 0.52476546298854465), ('c', 0.048283506458290146)])
      success: True
   parameters: array([  1.06371000e+00,   5.09155364e+04,   1.19022208e-06,
        -1.37443229e+00,  -2.74520718e-02])
   covariance: array([[  7.56396247e-01,  -1.47415535e-08,   9.27513342e-02,
          8.97227286e-03],
       [ -1.47415535e-08,   4.03514718e-15,  -1.57647791e-08,
          8.09694619e-10],
       [  9.27513342e-02,  -1.57647791e-08,   2.79008873e-01,
          1.00106837e-02],
       [  8.97227286e-03,   8.09694619e-10,   1.00106837e-02,
          2.35188405e-03]])
         ndof: 79
 vparam_names: array(['t0', 'x0', 'x1', 'c'], 
      dtype='|S2')
        chisq: 94.674280831892887
  param_names: array(['z', 't0', 'x0', 'x1', 'c'], 
      dtype='|S2')
      message: 'Minimization exited successfully.'
        ncall: 98,
 <sncosmo.models.Model at 0x108810ed0>)

In [14]:

    

iotest.mcmcOut


    

    
Out[14]:





(             vparam_names: array(['t0', 'x0', 'x1', 'c'], 
      dtype='|S2')
                   errors: OrderedDict([('t0', 0.84052291965928672), ('x0', 6.6097763429030765e-08), ('x1', 0.54376647284091717), ('c', 0.051716035384276604)])
                  samples: array([[  5.09172430e+04,   1.15286766e-06,  -5.73551900e-01,
          7.80925076e-02],
       [  5.09169532e+04,   1.15937672e-06,  -6.53440242e-01,
          6.60318272e-02],
       [  5.09169532e+04,   1.15937672e-06,  -6.53440242e-01,
          6.60318272e-02],
       ..., 
       [  5.09158980e+04,   1.22719862e-06,  -1.87470109e+00,
         -1.68128260e-02],
       [  5.09155473e+04,   1.24432644e-06,  -1.67843172e+00,
          4.68350001e-03],
       [  5.09155473e+04,   1.24432644e-06,  -1.67843172e+00,
          4.68350001e-03]])
               parameters: array([  1.06371000e+00,   5.09154497e+04,   1.19142568e-06,
        -1.26646862e+00,  -1.10709753e-02])
              param_names: array(['z', 't0', 'x0', 'x1', 'c'], 
      dtype='|S2')
 mean_acceptance_fraction: 0.59349999999999992
               covariance: array([[  7.06478778e-01,  -1.09831320e-08,   1.11635706e-01,
          9.94109082e-03],
       [ -1.09831320e-08,   4.36891433e-15,  -1.94720479e-08,
          9.84437375e-10],
       [  1.11635706e-01,  -1.94720479e-08,   2.95681977e-01,
          1.01745451e-02],
       [  9.94109082e-03,   9.84437375e-10,   1.01745451e-02,
          2.67454832e-03]]),
 <sncosmo.models.Model at 0x10953c890>)

In [15]:

    

iotest.nestOut


    

    
Out[15]:





(         logz: -60.441002229837324
         ndof: 79
        ncall: 6163
 vparam_names: array(['t0', 'x0', 'x1', 'c'], 
      dtype='|S2')
         logl: array([ -5.04238964e+04,  -3.96944503e+04,  -3.50742919e+04, ...,
        -4.74523891e+01,  -4.74136818e+01,  -4.74050143e+01])
   covariance: array([[  7.47158220e-01,  -1.22990009e-08,   6.04041064e-02,
          9.30440893e-03],
       [ -1.22990009e-08,   3.84070490e-15,  -1.49347510e-08,
          7.36404831e-10],
       [  6.04041064e-02,  -1.49347510e-08,   3.01669147e-01,
          1.04489323e-02],
       [  9.30440893e-03,   7.36404831e-10,   1.04489323e-02,
          2.46072809e-03]])
        niter: 1806
       errors: OrderedDict([('t0', 0.86438314449960862), ('x0', 6.197342094873836e-08), ('x1', 0.54924416000489951), ('c', 0.04960572640044595)])
   parameters: array([  1.06371000e+00,   5.09156021e+04,   1.18632748e-06,
        -1.26110938e+00,  -1.34699997e-02])
            h: 11.062512885719315
     logprior: array([ -4.61016602,  -4.62016602,  -4.63016602, ..., -22.66517019,
       -22.66517019, -22.66517019])
       bounds: {'x0': (0.0, 1.4249497311936984e-05), 'x1': (-3.0, 3.0), 'c': (-0.29999999999999999, 0.29999999999999999), 't0': (50774.025500000011, 51071.31319999999)}
   param_dict: OrderedDict([('t0', 50915.602091922599), ('x0', 1.1863274802416421e-06), ('x1', -1.2611093831758817), ('c', -0.013469999689044585), ('z', 1.0637099999999999)])
      weights: array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
         6.27371338e-05,   6.52131300e-05,   6.57808167e-05])
      samples: array([[  5.09150267e+04,   1.36766909e-05,   8.71338960e-01,
          2.42580968e-02],
       [  5.09108833e+04,   1.31580838e-05,  -1.04616311e+00,
         -1.34483275e-01],
       [  5.09755581e+04,   1.21386874e-05,   2.39309984e+00,
         -1.57536244e-01],
       ..., 
       [  5.09159287e+04,   1.18173883e-06,  -1.23987817e+00,
         -1.59102013e-02],
       [  5.09153033e+04,   1.17754694e-06,  -1.31047116e+00,
         -3.32303638e-02],
       [  5.09155558e+04,   1.19905724e-06,  -1.55125633e+00,
         -3.69537136e-02]])
         time: 398.26488494873047
      logzerr: 0.33260356110119016,
 <sncosmo.models.Model at 0x10dc859d0>)

In [13]:

    

figtest = test.plotLC()


    

In [17]:

    

figiotest = iotest.plotLC()


    

In [21]:

    

test.reset_all()


    

In [22]:

    

test.rerunFits()


    

    




mcmc_lc() output reran






    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-22-5f124b369968> in <module>()
----> 1 test.rerunFits()

/Users/lisaleemcb/sncosmo_lc_analysis/lc.py in rerunFits(self)
    155         if not self.mcmcRes:
    156             print "mcmc_lc() output reran"
--> 157             self.mcmcRes = self._mcmcOut[0]
    158             self.mcmcModel = self._mcmcOut[1]
    159 

TypeError: 'NoneType' object has no attribute '__getitem__'

In [18]:

    

test.plotCorner()


    

    
Out[18]:





(<matplotlib.figure.Figure at 0x109fa9d10>,
 <matplotlib.figure.Figure at 0x10c24fed0>)






    













    

In [19]:

    

iotest.plotCorner()


    

    
Out[19]:





(<matplotlib.figure.Figure at 0x10bc34810>,
 <matplotlib.figure.Figure at 0x112e7cf90>)