In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

In [2]:
import sncosmo

In [3]:
import glob

In [68]:
class TransientObject(object):
    
    def __init__(self, spectral_data, matched_spectra, days, data_fnames=None, mangled_fnames=None):
        """
        Constructor
        
        Parameters
        ----------
        spectral_data: list of 2-D `np.ndarray` of floats
            list of all the spectral data available corresponding to the object.
            Each `np.ndarray` has columns of wavelength and specific flux in units
            of Ang and arbitrary units respectively
        matched_spectra: list of 2-D `np.ndarray` of floats 
            list of all of the 'mangled' spectral data matched to the photometry of
            the transient object, in the form of 2-D `np.ndarray`. The list of arrays
            is indexed identically to the `spectral data`. Each `np.ndarray` has columns
            wavelength in units of Ang, and specific flux ($f_\lambda$) in units of 
            ergs/cm^2/s/Ang. It is assumed that the e
        days: `np.ndarray` of floats
            days on which the spectral_data were taken in the same 
        
        Returns
        -------
        instance of TransientObject class
        
        
        .. note :: 1. the wavelength gridding of the `mangled` matched spectra should 
        be identical, and may differ from the wavelength gridding of the spectral_data
        of the object from which it was constructed. 2. The wavelengths in the matched
        spectra are assumed to be in increasing order.
        """
        
        self.sortOrder = np.argsort(days)
        self.data_fnames = data_fnames[self.sortOrder]
        self.mangled_fnames = mangled_fnames[self.sortOrder]
        self.spectral_data = np.array(spectral_data)[self.sortOrder]
        self.mangled_spectra = np.array(matched_spectra)[self.sortOrder]
        self.mjds = np.array(days)[self.sortOrder]
        
        self.validateWavelengths()
        
        
    @classmethod
    def fromDataDir(cls, spectralDataDir):
        """
        """
        import glob

        
        spectralDataFname = np.array(glob.glob(spectralDataDir + '/*.dat'))
        mangledDataFname = np.array(glob.glob(spectralDataDir + '/*mangled.txt'))
            
        spectral_data = map(lambda x: np.loadtxt(x, skiprows=2), spectralDataFname) 
        matched_spectra = map(lambda x: np.loadtxt(x, skiprows=2), mangledDataFname)
        days = map(cls.parse_mjds, mangledDataFname)
        
        
            
        return cls(spectral_data=spectral_data, matched_spectra=matched_spectra, days=days,
                   data_fnames=spectralDataFname, mangled_fnames=mangledDataFname) 
            
    @staticmethod    
    def parse_mjds(mangledSpectraFileName):
        """
        
        """
        s = mangledSpectraFileName.split('_mangled')[0]
        day = s.split('_')[-1]
        return float(day)
    
    def validateWavelengths(self, verbose=False):
        """
        Check that the wavelength grid of the `self.matched_spectra
        are identical to each other and that the wavelengths are
        in increasing order
        """
        firstwave = self.mangled_spectra[0][:, 0]
        
        for i, mangledSpectra in enumerate(self.mangled_spectra):    
            wave = mangledSpectra[:, 0]
            if len(wave) != len(firstwave):
                print "Different lengths for ", 0, i
            elif not np.allclose(wave, firstwave):
                print 'Different values for ', 0, i
            else:
                if verbose:
                    print '0 ', i , ' match!'

In [8]:
location = '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy'

In [66]:
t  = TransientObject.fromDataDir(location)

In [67]:
t.mjds[np.argsort(t.mjds)]


Out[67]:
array([  0.     ,   1.9861 ,   6.95134,   7.94439,   9.93049,  25.8193 ,
        58.5899 ,  67.5273 ,  88.3813 ])

In [21]:
t.validateWavelengths()


0  0  match!
0  1  match!
0  2  match!
0  3  match!
0  4  match!
0  5  match!
0  6  match!
0  7  match!
0  8  match!

In [ ]:


In [22]:
import sncosmo

In [60]:
t.mangled_fnames


Out[60]:
array([ '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD0_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD1.9861_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD6.95134_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD7.94439_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD9.93049_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD25.8193_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD58.5899_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD67.5273_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD88.3813_mangled.txt',
       '/Users/rbiswas/soft/core_collapse_templates/data/2007uy/cfa_2007uy/cfa_2007uy_MJD153.923_mangled.txt'], 
      dtype='|S100')

In [29]:
fluxarray = np.zeros(shape=(len(t.mjds), len(t.mangled_spectra[0])))
for i, fluxdata in enumerate(t.mangled_spectra):
    fluxarray[i]= fluxdata[:, 1]

In [31]:
np.shape(fluxarray)


Out[31]:
(9, 2683)

In [43]:
for i in range(len(fluxarray)):
    plt.plot(t.mangled_spectra[0][:, 0], fluxarray[i])



In [38]:
len(fluxarray[0])


Out[38]:
2683

In [40]:
len(t.mangled_spectra[0])


Out[40]:
2683

In [49]:
np.all(np.diff(t.mangled_spectra[0][:, 0]) > 0)


Out[49]:
True

In [44]:
t.mjds - t.mjds[0]


Out[44]:
array([  0.     ,   1.9861 ,   6.95134,   7.94439,   9.93049,  25.8193 ,
        58.5899 ,  67.5273 ,  88.3813 ])

In [59]:
ts = sncosmo.TimeSeriesSource(phase=t.mjds - peak, wave=t.mangled_spectra[0][:, 0], flux=fluxarray)

In [60]:
peak = ts.peakphase('bessellB')

In [61]:
peak


Out[61]:
-1.4702144268534133e-14

In [62]:
model = sncosmo.Model(source=ts)

In [63]:
model.set(z=0.07)

In [64]:
print model


<Model at 0x10937e5d0>
source:
  class      : TimeSeriesSource
  name       : None
  version    : None
  phases     : [-9.55778, .., 78.8235] days
  wavelengths: [3454.82, .., 7369.95] Angstroms
parameters:
  z         = 0.070000000000000007
  t0        = 0.0
  amplitude = 1.0

In [ ]: