Python JLA_library.build_dictionary方法代码示例

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_date_of_max(options):
    import numpy
    from astropy.table import Table
    import JLA_library as JLA

    params=JLA.build_dictionary(options.config)

    # -----------  Read in the configuration file ------------


    lightCurveFits=JLA.get_full_path(params['lightCurveFits'])
    lightCurves=JLA.get_full_path(params['lightCurves'])
    adjlightCurves=JLA.get_full_path(params['adjLightCurves'])


    # ---------  Read in the list of SNe ---------------------
    SNe = Table.read(lightCurveFits, format='fits')

    nSNe=len(SNe)
    print 'There are %d SNe' % (nSNe)

    # -----------   The lightcurve fitting -------------------

    J=[]

    for SN in SNe:
        SNfile='lc-'+SN['name']+'.list'
        #print 'Examining %s' % SN['name']
        inputFile=lightCurves+SNfile
        outputFile=adjlightCurves+SNfile
        # If needed refit the lightcurve and insert the date of maximum into the input file
        JLA.insertDateOfMax(SN['name'].strip(),inputFile,outputFile,options.force)

    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_model(options):

    import numpy
    import astropy.io.fits as fits
    import JLA_library as JLA
    from astropy.table import Table
    from astropy.cosmology import FlatwCDM
    from scipy.interpolate import interp1d


    # -----------  Read in the configuration file ------------
    params=JLA.build_dictionary(options.config)

    # -----------  Read in the SN ordering ------------------------
    SNeList = numpy.genfromtxt(options.SNlist,
                               usecols=(0, 2),
                               dtype='S30,S200',
                               names=['id', 'lc'])
    nSNe = len(SNeList)

    for i, SN in enumerate(SNeList):
        SNeList['id'][i] = SNeList['id'][i].replace('lc-', '').replace('.list', '').replace('_smp', '')

    lcfile = JLA.get_full_path(params[options.lcfits])
    SNe = Table.read(lcfile, format='fits')

    print 'There are %d SNe' % (nSNe)

    indices = JLA.reindex_SNe(SNeList['id'], SNe)
    SNe = SNe[indices]

    redshift = SNe['zcmb']
    replace=(redshift < 0)

    # For SNe that do not have the CMB redshift
    redshift[replace]=SNe[replace]['zhel']
    print len(redshift)

    if options.raw:
        # Data from the bottom left hand figure of Mosher et al. 2014.
        # This is option ii) that is descibed above
        offsets=Table.read(JLA.get_full_path(params['modelOffset']),format='ascii.csv')
        Delta_M=interp1d(offsets['z'], offsets['offset'], kind='linear',bounds_error=False,fill_value='extrapolate')(redshift)
    else:
        Om_0=0.303 # JLA value in the wCDM model
        cosmo1 = FlatwCDM(name='SNLS3+WMAP7', H0=70.0, Om0=Om_0, w0=-1.0)
        cosmo2 = FlatwCDM(name='SNLS3+WMAP7', H0=70.0, Om0=Om_0, w0=-1.024)
        Delta_M=5*numpy.log10(cosmo1.luminosity_distance(redshift)/cosmo2.luminosity_distance(redshift))
    
    # Build the covariance matrix. Note that only magnitudes are affected
    Zero=numpy.zeros(nSNe)
    H=numpy.concatenate((Delta_M,Zero,Zero)).reshape(3,nSNe).ravel(order='F')
    C_model=numpy.matrix(H).T * numpy.matrix(H)

    date = JLA.get_date()
    fits.writeto('C_model_%s.fits' % (date),numpy.array(C_model),clobber=True) 

    return None

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_Cstat(options):
    """Python program to compute C_stat
    """

    import numpy
    import astropy.io.fits as fits
    from astropy.table import Table
    import JLA_library as JLA

    # -----------  Read in the configuration file ------------

    params=JLA.build_dictionary(options.config)

    # -----------  Read in the SN ordering ------------------------
    SNeList = numpy.genfromtxt(options.SNlist,
                               usecols=(0, 2),
                               dtype='S30,S200',
                               names=['id', 'lc'])
    nSNe = len(SNeList)

    for i, SN in enumerate(SNeList):
        SNeList['id'][i] = SNeList['id'][i].replace('lc-', '').replace('.list', '')

    lcfile = JLA.get_full_path(params[options.lcfits])
    SNe = Table.read(lcfile, format='fits')


    # -----------  Read in the data --------------------------

    print 'There are %d SNe in the sample' % (nSNe)

    indices = JLA.reindex_SNe(SNeList['id'], SNe)
    SNe=SNe[indices]

    C_stat=numpy.zeros(9*nSNe*nSNe).reshape(3*nSNe,3*nSNe)

    for i,SN in enumerate(SNe):
        cov=numpy.zeros(9).reshape(3,3)
        cov[0,0]=SN['dmb']**2.
        cov[1,1]=SN['dx1']**2.
        cov[2,2]=SN['dcolor']**2.
        cov[0,1]=SN['cov_m_s']
        cov[0,2]=SN['cov_m_c']
        cov[1,2]=SN['cov_s_c']
        # symmetrise
        cov=cov+cov.T-numpy.diag(cov.diagonal())
        C_stat[i*3:i*3+3,i*3:i*3+3]=cov

    # -----------  Read in the base matrix computed using salt2_stat.cc ------------

    if options.base!=None:
        C_stat+=fits.getdata(options.base)

    date = JLA.get_date()
    fits.writeto('C_stat_%s.fits' % date,C_stat,clobber=True) 

    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_dust(options):
    """Python program to compute C_dust
    """

    import numpy
    import astropy.io.fits as fits
    import os
    import JLA_library as JLA

    # ---------- Read in the SNe list -------------------------

    SNelist = numpy.genfromtxt(options.SNlist,
                               usecols=(0, 2),
                               dtype='S30,S110',
                               names=['id', 'lc'])

    for i, SN in enumerate(SNelist):
        SNelist['id'][i] = SNelist['id'][i].replace('lc-','').replace('.list','')

    # -----------  Read in the configuration file ------------

    params=JLA.build_dictionary(options.config)
    try:
        salt_path = JLA.get_full_path(params['defsaltModel'])
    except KeyError:
        salt_path = ''
        
    # -----------   The lightcurve fitting -------------------

    # Compute the offset between the nominal value of the extinciton 
    # and the adjusted value
    # We first compute the difference in light curve fit parameters for E(B-V) * (1+offset)
    offset = 0.1

    j = []

    for SN in SNelist:
        inputFile = SN['lc']
        print 'Fitting %s ' % (SN['lc'])
        workArea = JLA.get_full_path(options.workArea)
        dm, dx1, dc = JLA.compute_extinction_offset(SN['id'], inputFile, offset, workArea, salt_path)
        j.extend([dm, dx1, dc])
    
    # But we want to compute the impact of an offset that is twice as large, hence the factor of 4 in the expression
    # 2017/10/13
    # But we want to compute the impact of an offset that is half as large, hence the factor of 4 in the denominator
    # cdust = numpy.matrix(j).T * numpy.matrix(j) * 4.0
    cdust = numpy.matrix(j).T * numpy.matrix(j) / 4.0

    date = JLA.get_date()

    fits.writeto('C_dust_%s.fits' % date, cdust, clobber=True) 

    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_model(options):

    import numpy
    import astropy.io.fits as fits
    import JLA_library as JLA
    from astropy.table import Table
    from astropy.cosmology import FlatwCDM



    # -----------  Read in the configuration file ------------

    params=JLA.build_dictionary(options.config)

    # -----------  Read in the SN ordering ------------------------
    SNeList = numpy.genfromtxt(options.SNlist,
                               usecols=(0, 2),
                               dtype='S30,S200',
                               names=['id', 'lc'])
    nSNe = len(SNeList)

    for i, SN in enumerate(SNeList):
        SNeList['id'][i] = SNeList['id'][i].replace('lc-', '').replace('.list', '')

    lcfile = JLA.get_full_path(params[options.lcfits])
    SNe = Table.read(lcfile, format='fits')

    print 'There are %d SNe' % (nSNe)

    #z=numpy.array([])
    #offset=numpy.array([])
    Om_0=0.303 # JLA value in the wCDM model

    cosmo1 = FlatwCDM(name='SNLS3+WMAP7', H0=70.0, Om0=Om_0, w0=-1.0)
    cosmo2 = FlatwCDM(name='SNLS3+WMAP7', H0=70.0, Om0=Om_0, w0=-1.024)
    
    # For the JLA SNe
    redshift = SNe['zcmb']
    replace=(redshift < 0)
    # For the non JLA SNe
    redshift[replace]=SNe[replace]['zhel']

    Delta_M=5*numpy.log10(cosmo1.luminosity_distance(redshift)/cosmo2.luminosity_distance(redshift))

    # Build the covariance matrix. Note that only magnitudes are affected
    Zero=numpy.zeros(nSNe)
    H=numpy.concatenate((Delta_M,Zero,Zero)).reshape(3,nSNe).ravel(order='F')
    C_model=numpy.matrix(H).T * numpy.matrix(H)

    date = JLA.get_date()
    fits.writeto('C_model_%s.fits' % (date),numpy.array(C_model),clobber=True) 

    return None

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def convert_lightcurves(options):

    # Read in the configuration file
    # The configuraiton file contains the location of various files
    params=JLA.build_dictionary(options.config)

    # Read in the extra variance
    # This depends on the photometric method. It is lower for SMP 
    extraVariance=get_extra_variance(JLA.get_full_path(params['extraVariance']),options)

    # Read in the extinction values
    # A temporary fix as the lightcurves do not currently have it
    # Still needed
    extinction=get_extinction(JLA.get_full_path(params['extinction']),options)

    snanaDir=JLA.get_full_path(params['snanaLightCurves'])
    saltDir=JLA.get_full_path(params['adjLightCurves'])

    try:
        os.mkdir(saltDir)
    except:
        pass
        
    saltDir=saltDir+'DES/'

    try:
        os.mkdir(saltDir)
    except:
        pass

    for lightcurve in os.listdir(snanaDir):
        if '.dat' in lightcurve:
            # Read in the snana file
            lc=snanaLightCurve(snanaDir+lightcurve)
            lightCurveFile=saltDir+lightcurve.replace('des_real','lc-DES').replace('.dat','.list')
            if lc.parameters['TYPE'].split()[0] in ['1','101']:   # Is a SN Ia or a SN Ia?
                print lightcurve, lightCurveFile
                lc.clean()                                # Remove bad photometry
                lc.addNoise(extraVariance)                # Add additional variance to the lightcurve points
                # It is not clear if we need to compute a rough date of max before doing the more precise fit
                lc.estimateDateOfMax(options)             # Sets an approximate date of max for the light curve fitting done below.
                # Apply cuts
                # lc.applySNCuts()
                # lc.applySamplingCuts()
                lc.write(lightCurveFile,options.format)   # Write out the resutlt
                lc.fitDateOfMax(lightCurveFile,params)    # Get a more precise estimate of the data of peak brightness
                lc.updateExtinction(lightCurveFile,extinction) # Temporary code
    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_date_of_max(options):
    import numpy
    from astropy.table import Table
    import JLA_library as JLA

    params=JLA.build_dictionary(options.config)

    # ----------- Correction factor for extinction -----------
    # See ApJ 737 103
    extinctionFactor=0.86

    # -----------  Read in the configuration file ------------

    lightCurveFits=JLA.get_full_path(params['lightCurveFits'])
    lightCurves=JLA.get_full_path(params['lightCurves'])
    adjlightCurves=JLA.get_full_path(params['adjLightCurves'])

    # ---------  Read in the list of SNe ---------------------
    # One can either use an ASCII file with the SN list or a fits file
    if options.SNlist == None:
        SNe = Table.read(lightCurveFits, format='fits')
    else:
    # We use the ascii file, which gives the full path name
        SNe = Table.read(options.SNlist, format='ascii',names=['name','type','lc'],data_start=0)

    nSNe=len(SNe)
    print 'There are %d SNe' % (nSNe)

    # -----------   The lightcurve fitting -------------------

    for SN in SNe:
        if options.SNlist == None:
            SNfile='lc-'+SN['name']+'.list'
            inputFile=lightCurves+SNfile
            outputFile=adjlightCurves+SNfile
        else:
            inputFile=SN['lc']
            outputFile=SN['lc'].replace(lightCurves,adjlightCurves)

        print 'Examining %s' % SN['name']
        # If needed refit the lightcurve and insert the date of maximum into the input file
        JLA.insertDateOfMax(SN['name'].strip(),inputFile,outputFile,options.force,params)
        # Shouldn't we adjust the extinction first
        if options.adjustExtinction:
            adjustExtinction(outputFile,extinctionFactor)
    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_ZP(options):

    import JLA_library as JLA
    import numpy as np

    params=JLA.build_dictionary(options.config)

    # Read in the standard star

    standard=JLA.spectrum(JLA.get_full_path(params['magSys'])+options.standard)

    # Read in the filter

    filt=JLA.filterCurve(JLA.get_full_path(params['filterDir'])+options.filter)

    # Compute the ZP
    if options.system=='AB':
        print '%s in %s %s %5.3f' % (options.standard,options.filter,options.system,filt.AB(standard))
    else:
        pass
#        print '%s in %s %s %5.3f' % (options.standard,options.filter,options.system,filt.Vega(standard))


    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def create_Models(options):
    import os

    params=JLA.build_dictionary(options.config)

    try:
        os.mkdir(options.output)
    except:
        print "Directory %s already exists" % (options.output)

    # Read in the SALT models that will be kept
    SALTmodels=Table.read(options.modelList,format='ascii',names=['ID','Description'],data_start=0)

    # Read in the models for which the magnitude will be adjusted
    try:
        magOffsets=Table.read(options.magOffsetList,format='ascii',names=['Model','Filter','Offset','MagSys'],data_start=1,delimiter='\t',comment='#')
    except:
        magOffsets=[]

    modelList=[]
    for model in os.listdir(JLA.get_full_path(options.base)):
        if model in SALTmodels['ID']:
            print "Copying across %s" % model
            modelList.append(model)
            shutil.copytree(options.base+'/'+model,options.output+'/'+model)
            # Copy salt2 directory to salt2-4
            shutil.copytree(options.output+'/'+model+'/snfit_data/salt2',options.output+'/'+model+'/snfit_data/salt2-4')
            # Update fitmodel.card
            shutil.copy(JLA.get_full_path(params['fitmodel']),options.output+'/'+model+'/snfit_data/fitmodel.card')
            # Add the DECam instrument files
            shutil.copytree(JLA.get_full_path(params['DES_instrument']),options.output+'/'+model+'/snfit_data/Instruments/DECam')
            # Update the Keplercam instrument files
            shutil.rmtree(options.output+'/'+model+'/snfit_data/Instruments/Keplercam')
            shutil.copytree(JLA.get_full_path(params['CfA_instrument']),options.output+'/'+model+'/snfit_data/Instruments/Keplercam')
            ## Serious bug - filter and ZP offsets are lost here!
            # Add DES magnitude system
            shutil.copy(JLA.get_full_path(params['DES_magsys']),options.output+'/'+model+'/snfit_data/MagSys/')
            # Update the CfA magnitude system
            shutil.copy(JLA.get_full_path(params['CfA_magsys']),options.output+'/'+model+'/snfit_data/MagSys/')
            # This is not needed for CSP as the instrument files and magnitude system have not chnaged since JLA
        else:
            print "Excluding %s" % model


    print 'We start with %d models from JLA' % (len(modelList))

    # ---------  Add new models --------------

    newModels=Table.read(options.add,format='ascii', comment='#')
    for model in newModels:
        # Copy accross the base model
        shutil.copytree(JLA.get_full_path(model['baseModel']),options.output+'/'+model['modelNumber'])
        print 'Creating %s' % (model['modelNumber'])

        # Copy salt2 directory to salt2-4
        shutil.copytree(options.output+'/'+model['modelNumber']+'/snfit_data/salt2',options.output+'/'+model['modelNumber']+'/snfit_data/salt2-4')

        # Remove the old base instrument, if it exists and replace it with a new one
        try:
            shutil.rmtree(options.output+'/'+model['modelNumber']+'/snfit_data/'+model['fitmodel'])
        except:
            pass

        shutil.copytree(JLA.get_full_path(model['baseInstrument']+model['fitmodel']),options.output+'/'+model['modelNumber']+'/snfit_data/'+model['fitmodel'])

        # Remove the old MagSys directory and replace it with the new one
        shutil.rmtree(options.output+'/'+model['modelNumber']+'/snfit_data/MagSys')
        shutil.copytree(JLA.get_full_path(model['baseInstrument'])+'MagSys',options.output+'/'+model['modelNumber']+'/snfit_data/MagSys')
        
        # Replace the fitmodel.card it with the new one
        shutil.copy(JLA.get_full_path(model['baseInstrument'])+'fitmodel.card',options.output+'/'+model['modelNumber']+'/snfit_data/fitmodel.card')

        # Modify filter curve and ZP
        if model['Type']=='filt':
            offsetFilter(options.output+'/'+model['modelNumber']+'/snfit_data/'+model['fitmodel']+'/'+model['Filter'],model['Instrument'])
        else:
            offsetZP(options.output+'/'+model['modelNumber']+'/snfit_data/MagSys/'+model['MagSys'],model['ShortName'],model['Instrument'],model['fitmodel'])

        # We now update the list of instruments in the newly created surfaces
        # We should try to generalise this, as this will become very complex as more instruments are added.
        if model['Instrument']=='DECAM':
            # Update just the Keplercam instrument files
            updateKeplercam(options,params,model)
        elif model['Instrument']=='KEPLERCAM':
            # Update just the DES instrument files
            updateDES(options,params,model)
        else: # The case for swope ...
            # Update both the Keplercam and DES fles
            updateDES(options,params,model)
            updateKeplercam(options,params,model)


        modelList.append(model['modelNumber'])

    # ---- Update magnitude ZPs -----
    #for model in magOffsets:
    #    if numpy.abs(model['Offset']) > 0:
    #        magSys=options.output+'/'+model['Model']+'/snfit_data/MagSys/'+ model['MagSys']
    #        offsetZP2(magSys,model['Offset'],model['Filter'])

#.........这里部分代码省略.........

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def add_covar_matrices(options):
    """
    Python program that adds the individual covariance matrices into a single matrix
    """

    import time
    import numpy
    import astropy.io.fits as fits
    import JLA_library as JLA

    params = JLA.build_dictionary(options.config)

    # Read in the terms that account for uncertainties in perculiar velocities, 
    # instrinsic dispersion and, lensing

    # Read in the covariance matrices
    matrices = []

    systematic_terms = ['bias', 'cal', 'host', 'dust', 'model', 'nonia', 'pecvel', 'stat']

    covmatrices = {'bias':params['bias'],
                   'cal':params['cal'],
                   'host':params['host'],
                   'dust':params['dust'],
                   'model':params['model'],
                   'nonia':params['nonia'],
                   'pecvel':params['pecvel'],
                   'stat':params['stat']}

    for term in systematic_terms:
        matrices.append(fits.getdata(JLA.get_full_path(covmatrices[term]), 0))

    # Add the matrices
    size = matrices[0].shape
    add = numpy.zeros(size[0]**2.).reshape(size[0], size[0])
    for matrix in matrices:
        add += matrix

    # Write out this matrix. This is C_eta in qe. 13 of B14

    date=JLA.get_date()

    fits.writeto('C_eta_%s.fits' % (date), add, clobber=True)

    # Compute A

    nSNe = size[0]/3

    jla_results = {'Om':0.303, 'w':-1.027, 'alpha':0.141, 'beta':3.102}

    arr = numpy.zeros(nSNe*3*nSNe).reshape(nSNe, 3*nSNe)

    for i in range(nSNe):
        arr[i, 3*i] = 1.0
        arr[i, 3*i+1] = jla_results['alpha']
        arr[i, 3*i+2] = -jla_results['beta']

    cov = numpy.matrix(arr) * numpy.matrix(add) * numpy.matrix(arr).T

    # Add the diagonal terms

    sigma = numpy.genfromtxt(JLA.get_full_path(params['diag']),
                             comments='#',
                             usecols=(0, 1, 2),
                             dtype='f8,f8,f8',
                             names=['sigma_coh', 'sigma_lens', 'sigma_pecvel'])

    for i in range(nSNe):
        cov[i, i] += sigma['sigma_coh'][i]**2 + \
        sigma['sigma_lens'][i]**2 + \
        sigma['sigma_pecvel'][i]**2

    fits.writeto('C_total_%s.fits' % (date), cov, clobber=True)

    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_Ccal(options):
    """Python program to compute Ccal
    """

    import numpy
    import astropy.io.fits as fits
    from astropy.table import Table
    import multiprocessing as mp
    import matplotlib.pyplot as plt

    # -----------  Read in the configuration file ------------

    params=JLA.build_dictionary(options.config)
    try:
        salt_prefix = params['saltPrefix']
    except KeyError:
        salt_prefix = ''

    # ---------- Read in the SNe list -------------------------

    SNeList = Table(numpy.genfromtxt(options.SNlist,
                                     usecols=(0, 2),
                                     dtype='S30,S100',
                                     names=['id', 'lc']))


    for i,SN in enumerate(SNeList):
        SNeList['id'][i]=SNeList['id'][i].replace('lc-', '').replace('.list', '').replace('_smp', '')

    # ----------  Read in the SN light curve fits ------------
    # This is used to get the SN redshifts which are used in smoothing the Jacbian

    lcfile = JLA.get_full_path(params[options.lcfits])
    SNe = Table.read(lcfile, format='fits')

    # Make sure that the order is correct
    indices = JLA.reindex_SNe(SNeList['id'], SNe)
    SNe = SNe[indices]
    if len(indices) != len(SNeList['id']):
        print "We are missing SNe"
        exit()

    # -----------  Set up the structures to handle the different salt models -------
    # The first model is the unperturbed salt model
    SALTpath=JLA.get_full_path(params['saltPath'])

    SALTmodels=JLA.SALTmodels(SALTpath+'/saltModels.list')
    nSALTmodels=len(SALTmodels)-1
    print SALTmodels, nSALTmodels

    nSNe=len(SNeList)
    print 'There are %d SNe in the sample' % (nSNe)
    print 'There are %d SALT models' % (nSALTmodels)

    # Add a survey column, which we use with the smoothing, and the redshift
    SNeList['survey'] = numpy.zeros(nSNe,'a10')
    SNeList['z'] = SNe['zhel']

    # Identify the SNLS, SDSS, HST and low-z SNe. We use this when smoothing the Jacobian
    # There is rather inelegant 
    # We still need to allow for Vanina's naming convention when doing this for the photometric sample
    for i,SN in enumerate(SNeList):
        if SN['id'][0:4]=='SDSS':
            SNeList['survey'][i]='SDSS'
        elif SN['id'][2:4] in ['D1','D2','D3','D4']:
            SNeList['survey'][i]='SNLS'
        elif SN['id'][0:3] in ['DES']:
            SNeList['survey'][i]='DES'
        elif SN['id'][0:2]=='sn':
            SNeList['survey'][i]='nearby'
        else:
            SNeList['survey'][i]='high-z'

    # -----------   Read in the calibration matrix -----------------
    Cal=fits.getdata(JLA.get_full_path(params['C_kappa']))

    # Multiply the ZP submatrix by 100^2, and the two ZP-offset submatrices by 100,
    # because the magnitude offsets are 0.01 mag and the units of the covariance matrix are mag
    size=Cal.shape[0] / 2
    Cal[0:size,0:size]=Cal[0:size,0:size]*10000.
    Cal[0:size,size:]*=Cal[0:size,size:]*100.
    Cal[size:,0:size]=Cal[size:,0:size]*100.


    # ------------- Create an area to work in -----------------------
    workArea = JLA.get_full_path(options.workArea)
    try:
        os.mkdir(workArea)
    except:
        pass

    # -----------   The lightcurve fitting --------------------------

    firstSN=True
    
    log=open('log.txt','w')

    for i,SN in enumerate(SNeList):
        J=[]
        try:
#.........这里部分代码省略.........

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_bias(options):

    import numpy
    import astropy.io.fits as fits
    import JLA_library as JLA
    from astropy.table import Table
    from astropy.cosmology import FlatwCDM
    from  scipy.optimize import leastsq
    import matplotlib.pyplot as plt
    from scipy.stats import t


    # -----------  Read in the configuration file ------------
    params=JLA.build_dictionary(options.config)

    # -----------  Read in the SN ordering ------------------------
    SNeList = Table(numpy.genfromtxt(options.SNlist,
                               usecols=(0, 2),
                               dtype='S30,S200',
                               names=['id', 'lc']))
    nSNe = len(SNeList)

    for i, SN in enumerate(SNeList):
        SNeList['id'][i] = SNeList['id'][i].replace('lc-', '').replace('.list', '').replace('_smp','')

    lcfile = JLA.get_full_path(params[options.lcfits])
    SNe = Table.read(lcfile, format='fits')
    print 'There are %d SNe' % (nSNe)

    indices = JLA.reindex_SNe(SNeList['id'], SNe)
    SNe=SNe[indices]

    # Add a column that records the error in the bias correction
    SNe['e_bias'] = numpy.zeros(nSNe,'f8')

    # Read in the bias correction (see, for example, Fig.5 in B14)
    # Fit a polynomial to the data
    # Determine the uncertainties

    bias = numpy.genfromtxt(JLA.get_full_path(params['biasPolynomial']),
                                  skip_header=4,
                                  usecols=(0, 1, 2, 3),
                                  dtype='S10,f8,f8,f8',
                                  names=['sample', 'redshift', 'bias', 'e_bias'])

    if options.plot:
        fig=plt.figure()
        ax=fig.add_subplot(111)
        colour={'nearby':'b','SNLS':'r','SDSS':'g','DES':'k'}

    for sample in numpy.unique(bias['sample']):
        selection=(bias['sample']==sample)
        guess=[0,0,0]

        print bias[selection]
        plsq=leastsq(residuals, guess, args=(bias[selection]['bias'],
                                             bias[selection]['redshift'],
                                             bias[selection]['e_bias'],
                                             'poly'), full_output=1)

        if plsq[4] in [1,2,3,4]:
            print 'Solution for %s found' % (sample)

        if options.plot:
            ax.errorbar(bias[selection]['redshift'],
                    bias[selection]['bias'],
                    yerr=bias[selection]['e_bias'],
                    ecolor='k',
                    color=colour[sample],
                    fmt='o',
                    label=sample)
            z=numpy.arange(numpy.min(bias[selection]['redshift']),numpy.max(bias[selection]['redshift']),0.001)
            ax.plot(z,poly(z,plsq[0]),color=colour[sample])

        # For each SNe, determine the uncerainty in the correction. We use the approach descibed in
        # https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html
        
        # Compute the chi-sq.
        chisq=(((bias[selection]['bias']-poly(bias[selection]['redshift'],plsq[0]))/bias[selection]['e_bias'])**2.).sum()
        dof=selection.sum()-len(guess)
        print "Reduced chi-square value for sample %s is %5.2e" % (sample, chisq / dof)

        alpha=0.315 # Confidence interval is 100 * (1-alpha)
        # Compute the upper alpha/2 value for the student t distribution with dof
        thresh=t.ppf((1-alpha/2.0), dof)
        
        if options.plot and sample!='nearby':
            # The following is only valid for polynomial fitting functions, and we do not compute it for the nearby sample
            upper_curve=[]
            lower_curve=[]
            for x in z:
                vect=numpy.matrix([1,x,x**2.])
                offset=thresh * numpy.sqrt(chisq / dof * (vect*numpy.matrix(plsq[1])*vect.T)[0,0])
                upper_curve.append(poly(x,plsq[0])+offset)
                lower_curve.append(poly(x,plsq[0])-offset)

            ax.plot(z,lower_curve,'--',color=colour[sample])
            ax.plot(z,upper_curve,'--',color=colour[sample])

        # Compute the error in the bias
#.........这里部分代码省略.........

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compareSALTsurfaces(surface):

    import matplotlib.pyplot as plt

    # -----------  Read in the configuration file ------------

    params=JLA.build_dictionary(options.config)
    
    # -----------  Read in the SALT models -------------------

    surface1=readSALTsurface(JLA.get_full_path(params['model1'])+'salt2-4/')
#    surface2=readSALTsurface(JLA.get_full_path(params['model2'])+'salt2-4/')
    surface2=readSALTsurface(JLA.get_full_path(params['model2']))
    

    # -----------  Plot the surfaces ----------------------
    fig1=plt.figure()
    for axes,x1 in enumerate([-2,0,2]):
        ax=fig1.add_subplot(3,1,axes+1)
        flux1=surface1.template_0[options.phase]['flux'] + x1 * surface1.template_1[options.phase]['flux']
        flux2=surface2.template_0[options.phase]['flux'] + x1 * surface2.template_1[options.phase]['flux']
        ax.plot(surface1.template_0[options.phase]['wave'],flux1)
        ax.plot(surface2.template_0[options.phase]['wave'],flux2)
        ax.text(7000,0.3,"C=0 x1=%2d" % x1)

    ax.set_xlabel("wavelength ($\AA$)")
        
    plt.savefig(options.config.replace(".config","_SED.png"))

       
    # -----------  Plot the colour laws ----------------------

    # See salt2extinction.cc
    # Note the extrapolation

#    /*
#    ========================================================
#    VERSION 1
#    ========================================================
#    if(l_B<=l<=l_R)
#    ext = exp( color * constant * ( alpha*l + params(0)*l^2 + params(1)*l^3 + ... ))
#    = exp( color * constant * P(l) )
#    alpha = 1-params(0)-params(1)-...
#    if(l>l_R)
#    ext = exp( color * constant * ( P(l_R) + P'(l_R)*(l-l_R) ) )
#    if(l<l_B)
#    ext = exp( color * constant * ( P(l_B) + P'(l_B)*(l-l_B) ) )
#  
#    ======================================================== 
#    */

    constant=0.4 * numpy.log(10)
    fig3=plt.figure()
    ax3=fig3.add_subplot(111)
    wave=surface1.template_0[options.phase]['wave']
    wave_min=2800.
    wave_max=7000.
    
    wave_min_reduced=reduced_lambda(wave_min)
    wave_max_reduced=reduced_lambda(wave_max)
    
    # See salt2extinction.h
    reduced_wave=reduced_lambda(wave)

    # Model 1
    alpha1=1.0
    
    # There are 4 co-efficients in the colour law
    for coeff in surface1.colour_law['coeff']:
        alpha1-=coeff

    p1=numpy.zeros(len(reduced_wave))
        
    # Compute derivatives for extrapolations
    p1_derivative_min=derivative(alpha1,surface1,wave_min_reduced)
    p1_derivative_max=derivative(alpha1,surface1,wave_max_reduced)

    # Compute colour law at the points of extrapolations
    p1_wave_min_reduced=colourLaw(alpha1,surface1,wave_min_reduced)
    p1_wave_max_reduced=colourLaw(alpha1,surface1,wave_max_reduced)

    for index,rl in enumerate(reduced_wave):
        if rl < wave_min_reduced:
            p1[index]=p1_wave_min_reduced+p1_derivative_min*(rl-wave_min_reduced)
        elif rl > wave_max_reduced:
            p1[index]=p1_wave_max_reduced+p1_derivative_max*(rl-wave_max_reduced)
        else:
            p1[index]=colourLaw(alpha1,surface1,rl)
    
    # Model 2
    alpha2=1.0

    for coeff in surface2.colour_law['coeff']:
        alpha2-=coeff

    p2=numpy.zeros(len(reduced_wave))
            
    # Compute derivatives for extrapolations
    p2_derivative_min=derivative(alpha2,surface2,wave_min_reduced)
    p2_derivative_max=derivative(alpha2,surface2,wave_max_reduced)
#.........这里部分代码省略.........

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_nonIa(options):
    """Pythom program to compute the systematic unsertainty related to
    the contamimation from Ibc SNe"""

    import numpy
    import astropy.io.fits as fits
    from astropy.table import Table, MaskedColumn, vstack
    import JLA_library as JLA

    # The program computes the covaraince for the spectroscopically confirmed SNe Ia only
    # The prgram assumes that the JLA SNe are first in any list
    # Taken from C11

    # Inputs are the rates of SNe Ia and Ibc, the most likely contaminant

    # Ia rate - Perett et al.
    # SN Ibc rate - proportional to the star formation rate - Hopkins and Beacom
    # SN Ib luminosity distribution. Li et al + bright SN Ibc Richardson

    # The bright Ibc population
    # d_bc = 0.25     # The offset in magnitude between the Ia and bright Ibc
    # s_bc = 0.25     # The magnitude spread
    # f_bright = 0.25 # The fraction of Ibc SN that are bright

    # Simulate the characteristics of the SNLS survey
    # Apply outlier rejection
    # All SNe that pass the cuts are included in the sample

    # One then has a mixture of SNe Ia and SNe Ibc
    # and the average magnitude at each redshift is biased. This
    # is called the raw bias. One multiplies the raw bias by the fraction of
    # objects classified as SNe Ia*

    # The results are presented in 7 redshift bins defined in table 14 of C11
    # We use these results to generate the matrix.
    # Only the SNLS SNe in the JLA sample are considered.
    # For the photometrically selected sample and other surveys, this will probably be different
    # JLA compute this for the SNLS sample only

    # We assume that the redshift in this table refers to the left hand edge of each bin

    # -----------  Read in the configuration file ------------

    params = JLA.build_dictionary(options.config)
    

    data=numpy.genfromtxt(JLA.get_full_path(params['classification']),comments="#",usecols=(0,1,2),dtype=['float','float','float'],names=['redshift','raw_bias','fraction'])
    z_bin=data['redshift']
    raw_bias=data['raw_bias']
    f_star=data['fraction']
    
    # The covaraiance between SNe Ia in the same redshift bin is fully correlated
    # Otherwise, it is uncorrelated

    # -----------  Read in the configuration file ------------

    params = JLA.build_dictionary(options.config)

    SNeList = numpy.genfromtxt(options.SNlist,
                               usecols=(0, 2),
                               dtype='S30,S200',
                               names=['id', 'lc'])

    for i, SN in enumerate(SNeList):
        SNeList['id'][i] = SNeList['id'][i].replace('lc-', '').replace('.list', '').replace('_smp','')

    lcfile = JLA.get_full_path(params[options.lcfits])
    SNe = Table.read(lcfile, format='fits')

    # Add a bin column and a column that specifies if the covariance needs to be computed
    SNe['bin'] = 0
    SNe['eval'] = False

    # make the order of data (in SNe) match SNeList
    indices = JLA.reindex_SNe(SNeList['id'], SNe)
    SNe = SNe[indices]

    nSNe = len(SNe)
    # Identify the SNLS SNe in the JLA sample
    # We use the source and the name to decide if we want to add corrections for non-Ia contamination
    # Identify the DESS SNe in the DES sample.
    for i, SN in enumerate(SNe):
        try:
            # If the source keyword exists
            if (SN['source'] == 'JLA' or SN['source'] == 'SNLS_spec') and SN['name'][2:4] in ['D1', 'D2', 'D3', 'D4']:
                SNe['eval'][i] = True
            elif (SN['source']== 'SNLS_photo') and (SN['name'][2:4] in ['D1', 'D2', 'D3', 'D4'] or (SN['name'][0:2] in ['D1', 'D2', 'D3', 'D4'])):
                SNe['eval'][i] = True
        except:
            # If the source keyword does not exist
            if SN['name'][0:3]=="DES":
                SNe['eval'][i] = True

    print list(SNe['eval']).count(True)
    # Work out which redshift bin each SNe belongs to
    # In numpy.digitize, the bin number starts at 1, so we subtract 1 -- need to check...
    SNe['bin'] = numpy.digitize(SNe['zhel'], z_bin)-1 

    # Build the covariance matrix
    C_nonIa = numpy.zeros(nSNe*3*nSNe*3).reshape(nSNe*3, nSNe*3)
#.........这里部分代码省略.........

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import build_dictionary [as 别名]
def compute_nonIa_fraction(options):

    import numpy
    import astropy.io.fits as fits
    from astropy.table import Table
    import JLA_library as JLA
    import matplotlib.pyplot as plt

    # -----------  Read in the configuration file ------------

    params = JLA.build_dictionary(options.config)

    # -----------  Read in the SNe ---------------------------

    SNeList = numpy.genfromtxt(options.SNlist,
                               usecols=(0, 2),
                               dtype='S30,S200',
                               names=['id', 'lc'])

    for i, SN in enumerate(SNeList):
        SNeList['id'][i] = SNeList['id'][i].replace('lc-', '').replace('.list', '')


    redshift=[]
    SNtype=[]
    for SN in SNeList:
        if 'DES' in SN['lc']:
            z,t=getVariable(SN['lc'],['Z_HELIO','SNTYPE'])
            redshift.append(float(z))
            SNtype.append(int(t))
            nSNe=len(redshift)
    types=numpy.zeros(nSNe,dtype=[('redshift','f4'),('SNtype','i4')])
    types['redshift']=redshift
    types['SNtype']=SNtype


    print "There are %d SNe" % (nSNe)



    bins=numpy.array([0.00,0.10,0.26,0.41,0.57,0.72,0.89,1.04])
    bias=numpy.array([0.00,0.015,0.024,0.024,0.024,0.023,0.026,0.025])
    selection=(types['SNtype']==1)
    SNeIa=numpy.histogram(types[selection]['redshift'],bins)
    SNeIa_total=numpy.histogram(types['redshift'],bins)
    print SNeIa, SNeIa_total

    if options.plot:
        fig=plt.figure()
        ax=fig.add_subplot(111)
        ax.hist(types[selection]['redshift'],bins=bins, histtype='step', color='b')
        ax.hist(types['redshift'],bins=bins, histtype='step', color='g')
        plt.show()

    if nSNe > SNeIa_total[0].sum():
        print 'Oops'
        exit()

    # Write out as an ascii table
    output=open(JLA.get_full_path(params['classification']),'w')
    output.write('# DES SN classification\n')
    output.write('# Written on %s\n' % (JLA.get_date()))
    output.write('# redshift\tRaw_bias\tfraction\n')
    output.write('%4.2f\t%4.3f\t%4.3f\n' % (0.0,0.0,0.0))
    for i in range(len(bins)-1):
        if SNeIa_total[0][i] > 0:
            output.write('%4.2f\t%4.3f\t%4.3f\n' % (bins[i+1],bias[i+1],1.0-1.0*SNeIa[0][i]/SNeIa_total[0][i]))

    output.close()
    

    return