Python JLA_library.get_full_path方法代码示例

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import get_full_path [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 get_full_path [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 get_full_path [as 别名]
def add_covar_matrices(covmatrices,diag):
    """
    Python program that adds the individual covariance matrices into a single matrix
    """

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

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

    # Read in the covariance matrices
    matrices = []
    for matrix in covmatrices:
        matrices.append(fits.getdata(JLA.get_full_path(covmatrices[matrix]), 0))
        # Test for NaNs and replace them with zero
        if numpy.isnan(matrices[-1]).any():
            print 'Found a NaN in %s ... replacing them with zero' % (covmatrices[matrix])
            print numpy.isnan(matrices[-1]).sum()
            matrices[-1][numpy.isnan(matrices[-1])]=0.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

    # 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(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

    return cov

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import get_full_path [as 别名]
def updateDES(options,params,model):
    try:
        shutil.rmtree(options.output+'/'+model['modelNumber']+'/snfit_data/Instruments/DECam')
    except:
        pass
    shutil.copytree(JLA.get_full_path(params['DES_instrument']),options.output+'/'+model['modelNumber']+'/snfit_data/Instruments/DECam')
    
    # Update the DES magnitude system
    shutil.copy(JLA.get_full_path(params['DES_magsys']),options.output+'/'+model['modelNumber']+'/snfit_data/MagSys/')

    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import get_full_path [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 get_full_path [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 get_full_path [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 get_full_path [as 别名]
def updateKeplercam(options,params,model):
    try:
        shutil.rmtree(options.output+'/'+model['modelNumber']+'/snfit_data/Instruments/Keplercam')
    except:
        pass
    shutil.copytree(JLA.get_full_path(params['KelplerCam_instrument']),options.output+'/'+model['modelNumber']+'/snfit_data/Instruments/Keplercam')
    
    return

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import get_full_path [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 get_full_path [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 get_full_path [as 别名]
def runSALT(SALTpath, SALTmodel, salt_prefix, inputFile, SN):
    import os
    
    # Set up the path to the SALT model and the name of the outputFile
    os.environ['SALTPATH']=SALTpath+SALTmodel['directory']+'/snfit_data/'
    outputFile=JLA.get_full_path(options.workArea)+'/'+SN+'/'+SN+'_'+SALTmodel['directory']+'.dat'
    if os.path.isfile(outputFile):
        pass
        #print "Skipping, fit with SALT model %s for %s already done" % (SALTmodel['directory'],os.path.split(inputFile)[1])
    else:
        # Otherwise, do the fit with the date of Max set to the value in the lightcurve file
        JLA.fitLC(inputFile, outputFile, salt_prefix, forceDayMax=True)
    return outputFile

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import get_full_path [as 别名]
def prop_unc(params,filt,spectrum=None):

    # Use the filterwheel to find the filename of the filter
    filterDir=params['DES_instrument']
    filterWheel=JLA.get_full_path(filterDir)+'/Filterwheel'
    filterNames=numpy.genfromtxt(filterWheel,comments='#',usecols=(0,1),dtype='S1,S30',
                                 names=['filterName','filename'])
    filter_filename=filterNames[filterNames['filterName']==filt['filter'][-1:]]['filename'][0]

    # Read in the filter curve
    filterCurve=JLA.filterCurve(JLA.get_full_path(filterDir)+'/'+filter_filename)

    # Set the magnitude of the filter offset
    offset=filt['wavelength']*10.  

    # We compute a number of integrals. First with the filtercurves as is, then with an offset added to the filter curve
    # i) The I0 integral 
    error_I0=2.5 * numpy.log10(filterCurve.I0(0.0)/filterCurve.I0(offset))
    # ii) The chromatic correction.
    # Assumed to be zero for now
    # If the standard filter transmission curves are shifted by 5nm, then all the filters will be out by that same amount
    # This may mean, that the offset is quite small
    #mean_wavelength=filterCurve.mean()
    #I10_std=filterCurve.I1(mean_wavelength,0.0) / filterCurve.I0(0.0)
    #I10_std_offset=filterCurve.I1(mean_wavelength,10.0) / filterCurve.I0(10.0)

    error_chromatic=0.0

    # iii) The error in the AB offset
    # We use the standard filter curve to compute the AB offset
    if spectrum==None:
        calspec=JLA.spectrum(fits.getdata(JLA.get_full_path(params['calspec']),1),'CALSPEC')
    else:
        calspec=JLA.spectrum(fits.getdata(JLA.get_full_path(spectrum),1),'CALSPEC')

    error_AB=filterCurve.AB(calspec)-filterCurve.AB(calspec,offset)
    return error_I0,error_chromatic,error_AB

# 需要导入模块: import JLA_library [as 别名]
# 或者: from JLA_library import get_full_path [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 get_full_path [as 别名]
 def fitDateOfMax(self,lightCurveFile,params):
     # A full salt2 fit
     outputFile=lightCurveFile.replace('.list','.res')
     os.environ['SALTPATH']=JLA.get_full_path(params['defsaltModel'])
     JLA.fitLC(lightCurveFile, outputFile, salt_prefix='')
     self.dateofMax,self.dateofMaxError=JLA.getDateOfMax(outputFile)
     # Remove the old date of max and insert the new one
     lc=open(lightCurveFile)
     lc_lines=lc.readlines()
     lc.close()
     lc=open(lightCurveFile,'w')
     lc.write('@DayMax %s %s\n' % (self.dateofMax,self.dateofMaxError))
     for line in lc_lines:
         if 'DayMax' in line:
             pass
         else:
             lc.write(line)
     lc.close()
     return

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

    eff=[]
    # Read in the filter curves
    filt=Table.read(JLA.get_full_path(options.filter),format='ascii.csv')
    for col in filt.colnames:
        if 'ccd' in col and 'amp' not in col:
            f=JLA.filterCurve(filt['wavelength'],filt[col])
            eff.append(f.eff())

    print 'Examined the transmission curves for %d CCDs' % (len(eff))
    print 'Mean effective wavelength is %6.1f' % (numpy.mean(eff))
    print 'Range of effective wavelength is %6.1f-%6.1f' % (numpy.min(eff),numpy.max(eff))
    print 'RMS effective wavelength %6.1f' % (numpy.std(eff))

    return