Python Sed.Sed類代碼示例

    def calcBasicColors(self, sedList, bandpassDict, makeCopy = False):

        """
        This will calculate a set of colors from a list of SED objects when there is no need to redshift
        the SEDs.

        @param [in] sedList is the set of spectral objects from the models SEDs provided by loaders in
        rgStar or rgGalaxy. NOTE: Since this uses photometryBase.manyMagCalc_list the SED objects
        will be changed.

        @param [in] bandpassDict is a BandpassDict class instance with the Bandpasses set to those
        for the magnitudes given for the catalog object

        @param [in] makeCopy indicates whether or not to operate on copies of the SED objects in sedList
        since this method will change the wavelength grid.

        @param [out] modelColors is the set of colors in the Bandpasses provided for the given sedList.
        """

        modelColors = []

        for specObj in sedList:
            if makeCopy==True:
                fileSED = Sed()
                fileSED.setSED(wavelen = specObj.wavelen, flambda = specObj.flambda)
                sEDMags = bandpassDict.magListForSed(fileSED)
            else:
                sEDMags = bandpassDict.magListForSed(specObj)
            colorInfo = []
            for filtNum in range(0, len(bandpassDict)-1):
                colorInfo.append(sEDMags[filtNum] - sEDMags[filtNum+1])
            modelColors.append(colorInfo)

        return modelColors

    def testAlternateBandpassesStars(self):
        """
        This will test our ability to do photometry using non-LSST bandpasses.

        It will first calculate the magnitudes using the getters in cartoonPhotometryStars.

        It will then load the alternate bandpass files 'by hand' and re-calculate the magnitudes
        and make sure that the magnitude values agree.  This is guarding against the possibility
        that some default value did not change and the code actually ended up loading the
        LSST bandpasses.
        """

        obs_metadata_pointed = ObservationMetaData(mjd=2013.23,
                                                   boundType='circle',
                                                   pointingRA=200.0, pointingDec=-30.0,
                                                   boundLength=1.0)

        test_cat = cartoonStars(self.star, obs_metadata=obs_metadata_pointed)

        with lsst.utils.tests.getTempFilePath('.txt') as catName:
            test_cat.write_catalog(catName)
            with open(catName, 'r') as input_file:
                lines = input_file.readlines()
                self.assertGreater(len(lines), 1)

        cartoonDir = os.path.join(getPackageDir('sims_photUtils'), 'tests', 'cartoonSedTestData')
        testBandPasses = {}
        keys = ['u', 'g', 'r', 'i', 'z']

        bplist = []

        for kk in keys:
            testBandPasses[kk] = Bandpass()
            testBandPasses[kk].readThroughput(os.path.join(cartoonDir, "test_bandpass_%s.dat" % kk))
            bplist.append(testBandPasses[kk])

        sedObj = Sed()
        phiArray, waveLenStep = sedObj.setupPhiArray(bplist)

        i = 0

        # since all of the SEDs in the cartoon database are the same, just test on the first
        # if we ever include more SEDs, this can be something like
        # for ss in test_cata.sedMasterList:
        ss = test_cat.sedMasterList[0]
        ss.resampleSED(wavelen_match = bplist[0].wavelen)
        ss.flambdaTofnu()
        mags = -2.5*np.log10(np.sum(phiArray*ss.fnu, axis=1)*waveLenStep) - ss.zp
        self.assertEqual(len(mags), len(test_cat.cartoonBandpassDict))
        self.assertGreater(len(mags), 0)
        for j in range(len(mags)):
            self.assertAlmostEqual(mags[j], test_cat.magnitudeMasterList[i][j], 4)

    def loadwdSEDs(self, subset = None):

        """
        By default will load all seds in wd directory. The user can also define a subset of
        what's in the directory and load only those SEDs instead. Will skip over extraneous
        files in sed folder.

        @param [in] subset is the list of the subset of files wanted if one doesn't want all files
        in the kurucz directory.

        @param [out] sedListH is the set of model SED spectra objects for Hydrogen WDs to be passed onto
        the matching routines.

        @param [out] sedListHE is the set of model SED spectra objects for Helium WDs to be passed onto
        the matching routines.
        """
        files = []

        if subset is None:
            for fileName in os.listdir(self.wdDir):
                files.append(fileName)
        else:
            for fileName in subset:
                files.append(fileName)

        numFiles = len(files)
        numOn = 0

        sedListH = []
        sedListHE = []

        for fileName in files:
            if numOn % 100 == 0:
                print 'Loading %i of %i: WD SEDs' % (numOn, numFiles)

            try:
                spec = Sed()
                spec.readSED_flambda(str(self.wdDir + '/' + fileName))
                spec.name = fileName
                if fileName.split("_")[1] == 'He':
                    sedListHE.append(spec)
                else:
                    sedListH.append(spec)

            except:
                continue

            numOn += 1

        return sedListH, sedListHE

    def testAlternateBandpassesStars(self):
        """
        This will test our ability to do photometry using non-LSST bandpasses.

        It will first calculate the magnitudes using the getters in cartoonPhotometryStars.

        It will then load the alternate bandpass files 'by hand' and re-calculate the magnitudes
        and make sure that the magnitude values agree.  This is guarding against the possibility
        that some default value did not change and the code actually ended up loading the
        LSST bandpasses.
        """

        obs_metadata_pointed=ObservationMetaData(mjd=2013.23,
                                                 boundType='circle',unrefractedRA=200.0,unrefractedDec=-30.0,
                                                 boundLength=1.0)

        test_cat=cartoonStars(self.star,obs_metadata=obs_metadata_pointed)
        test_cat.write_catalog("testStarsCartoon.txt")

        cartoonDir = os.getenv('SIMS_PHOTUTILS_DIR')+'/tests/cartoonSedTestData/'
        testBandPasses = {}
        keys = ['u','g','r','i','z']

        bplist = []

        for kk in keys:
            testBandPasses[kk] = Bandpass()
            testBandPasses[kk].readThroughput(os.path.join(cartoonDir,"test_bandpass_%s.dat" % kk))
            bplist.append(testBandPasses[kk])

        sedObj = Sed()
        phiArray, waveLenStep = sedObj.setupPhiArray(bplist)

        i = 0

        #since all of the SEDs in the cartoon database are the same, just test on the first
        #if we ever include more SEDs, this can be something like
        #for ss in test_cata.sedMasterList:
        #
        ss=test_cat.sedMasterList[0]
        ss.resampleSED(wavelen_match = bplist[0].wavelen)
        ss.flambdaTofnu()
        mags = -2.5*numpy.log10(numpy.sum(phiArray*ss.fnu, axis=1)*waveLenStep) - ss.zp
        self.assertTrue(len(mags)==len(test_cat.bandpassDict))
        self.assertTrue(len(mags)>0)
        for j in range(len(mags)):
            self.assertAlmostEqual(mags[j],test_cat.magnitudeMasterList[i][j],10)
        i += 1

        os.unlink("testStarsCartoon.txt")

    def testSignalToNoise(self):
        """
        Test that calcSNR_m5 and calcSNR_sed give similar results
        """
        defaults = LSSTdefaults()
        photParams = PhotometricParameters()
        totalDict, hardwareDict = BandpassDict.loadBandpassesFromFiles()

        skySED = Sed()
        skySED.readSED_flambda(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "darksky.dat"))

        m5 = []
        for filt in totalDict:
            m5.append(calcM5(skySED, totalDict[filt], hardwareDict[filt], photParams, seeing=defaults.seeing(filt)))

        sedDir = lsst.utils.getPackageDir("sims_sed_library")
        sedDir = os.path.join(sedDir, "starSED", "kurucz")
        fileNameList = os.listdir(sedDir)

        numpy.random.seed(42)
        offset = numpy.random.random_sample(len(fileNameList)) * 2.0

        for ix, name in enumerate(fileNameList):
            if ix > 100:
                break
            spectrum = Sed()
            spectrum.readSED_flambda(os.path.join(sedDir, name))
            ff = spectrum.calcFluxNorm(m5[2] - offset[ix], totalDict.values()[2])
            spectrum.multiplyFluxNorm(ff)
            magList = []
            controlList = []
            magList = []
            for filt in totalDict:
                controlList.append(
                    calcSNR_sed(
                        spectrum, totalDict[filt], skySED, hardwareDict[filt], photParams, defaults.seeing(filt)
                    )
                )

                magList.append(spectrum.calcMag(totalDict[filt]))

            testList, gammaList = calcSNR_m5(
                numpy.array(magList), numpy.array(totalDict.values()), numpy.array(m5), photParams
            )

            for tt, cc in zip(controlList, testList):
                msg = "%e != %e " % (tt, cc)
                self.assertTrue(numpy.abs(tt / cc - 1.0) < 0.001, msg=msg)

    def testSystematicUncertainty(self):
        """
        Test that systematic uncertainty is added correctly.
        """
        sigmaSys = 0.002
        m5 = [23.5, 24.3, 22.1, 20.0, 19.5, 21.7]
        photParams = PhotometricParameters(sigmaSys=sigmaSys)

        bandpassDict = BandpassDict.loadTotalBandpassesFromFiles()
        obs_metadata = ObservationMetaData(unrefractedRA=23.0, unrefractedDec=45.0, m5=m5, bandpassName=self.bandpasses)
        magnitudes = bandpassDict.magListForSed(self.starSED)

        skySeds = []

        for i in range(len(self.bandpasses)):
            skyDummy = Sed()
            skyDummy.readSED_flambda(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "darksky.dat"))
            normalizedSkyDummy = setM5(
                obs_metadata.m5[self.bandpasses[i]],
                skyDummy,
                self.totalBandpasses[i],
                self.hardwareBandpasses[i],
                seeing=LSSTdefaults().seeing(self.bandpasses[i]),
                photParams=PhotometricParameters(),
            )

            skySeds.append(normalizedSkyDummy)

        for i in range(len(self.bandpasses)):
            snr = calcSNR_sed(
                self.starSED,
                self.totalBandpasses[i],
                skySeds[i],
                self.hardwareBandpasses[i],
                seeing=LSSTdefaults().seeing(self.bandpasses[i]),
                photParams=PhotometricParameters(),
            )

            testSNR, gamma = calcSNR_m5(
                numpy.array([magnitudes[i]]),
                [self.totalBandpasses[i]],
                numpy.array([m5[i]]),
                photParams=PhotometricParameters(sigmaSys=0.0),
            )

            self.assertAlmostEqual(snr, testSNR[0], 10, msg="failed on calcSNR_m5 test %e != %e " % (snr, testSNR[0]))

            control = numpy.sqrt(numpy.power(magErrorFromSNR(testSNR), 2) + numpy.power(sigmaSys, 2))

    def __init__(self, catsim_cat, om10_cat='twinkles_lenses_v2.fits',
                 density_param=1.):
        """
        Parameters
        ----------
        catsim_cat: catsim catalog
            The results array from an instance catalog.
        om10_cat: optional, defaults to 'twinkles_tdc_rung4.fits
            fits file with OM10 catalog
        density_param: `np.float`, optioanl, defaults to 1.0
            the fraction of eligible agn objects that become lensed and should
            be between 0.0 and 1.0.

        Returns
        -------
        updated_catalog:
            A new results array with lens systems added.
        """

        twinklesDir = getPackageDir('Twinkles')
        om10_cat = os.path.join(twinklesDir, 'data', om10_cat)
        self.catalog = catsim_cat
        # ****** THIS ASSUMES THAT THE ENVIRONMENT VARIABLE OM10_DIR IS SET *******
        lensdb = om10.DB(catalog=om10_cat)
        self.lenscat = lensdb.lenses.copy()
        self.density_param = density_param
        self.bandpassDict = BandpassDict.loadTotalBandpassesFromFiles(bandpassNames=['i'])

        specFileStart = 'Burst'
        for key, val in sorted(iteritems(SpecMap.subdir_map)):
            if re.match(key, specFileStart):
                galSpecDir = str(val)
        galDir = str(getPackageDir('sims_sed_library') + '/' + galSpecDir + '/')
        self.LRG_name = 'Burst.25E09.1Z.spec'
        self.LRG = Sed()
        self.LRG.readSED_flambda(str(galDir + self.LRG_name))
        #return

        #Calculate imsimband magnitudes of source galaxies for matching
        agn_sed = Sed()
        agn_fname = str(getPackageDir('sims_sed_library') + '/agnSED/agn.spec.gz')
        agn_sed.readSED_flambda(agn_fname)
        src_iband = self.lenscat['MAGI_IN']
        self.src_mag_norm = []
        for src in src_iband:
            self.src_mag_norm.append(matchBase().calcMagNorm([src],
                                                             agn_sed,
                                                             self.bandpassDict))

    def __init__(self, catsim_cat, om10_cat='twinkles_tdc_rung4.fits', density_param = 1.):
        """
        Input:
        catsim_cat:
            The results array from an instance catalog.

        density_param:
            A float between 0. and 1.0 that determines the fraction of eligible agn objects that become lensed.

        Output:
        updated_catalog:
            A new results array with lens systems added.
        """


        self.catalog = catsim_cat
        # ****** THIS ASSUMES THAT THE ENVIRONMENT VARIABLE OM10_DIR IS SET *******
        lensdb = om10.DB(catalog=om10_cat)
        self.lenscat = lensdb.lenses.copy()
        self.density_param = density_param
        self.bandpassDict = BandpassDict.loadTotalBandpassesFromFiles(bandpassNames=['i'])

        specFileStart = 'Burst'
        for key, val in sorted(iteritems(SpecMap.subdir_map)):
            if re.match(key, specFileStart):
                galSpecDir = str(val)
        galDir = str(getPackageDir('sims_sed_library') + '/' + galSpecDir + '/')
        self.LRG_name = 'Burst.25E09.1Z.spec'
        self.LRG = Sed()
        self.LRG.readSED_flambda(str(galDir + self.LRG_name))

    def setUp(self):
        starName = os.path.join(lsst.utils.getPackageDir("sims_sed_library"), defaultSpecMap["km20_5750.fits_g40_5790"])
        self.starSED = Sed()
        self.starSED.readSED_flambda(starName)
        imsimband = Bandpass()
        imsimband.imsimBandpass()
        fNorm = self.starSED.calcFluxNorm(22.0, imsimband)
        self.starSED.multiplyFluxNorm(fNorm)

        self.totalBandpasses = []
        self.hardwareBandpasses = []

        componentList = ["detector.dat", "m1.dat", "m2.dat", "m3.dat", "lens1.dat", "lens2.dat", "lens3.dat"]
        hardwareComponents = []
        for c in componentList:
            hardwareComponents.append(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", c))

        self.bandpasses = ["u", "g", "r", "i", "z", "y"]
        for b in self.bandpasses:
            filterName = os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "filter_%s.dat" % b)
            components = hardwareComponents + [filterName]
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughputList(components)
            self.hardwareBandpasses.append(bandpassDummy)
            components = components + [os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "atmos.dat")]
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughputList(components)
            self.totalBandpasses.append(bandpassDummy)

    def setUp(self):
        starName = os.path.join(lsst.utils.getPackageDir('sims_sed_library'),defaultSpecMap['km20_5750.fits_g40_5790'])
        self.starSED = Sed()
        self.starSED.readSED_flambda(starName)
        imsimband = Bandpass()
        imsimband.imsimBandpass()
        fNorm = self.starSED.calcFluxNorm(22.0, imsimband)
        self.starSED.multiplyFluxNorm(fNorm)

        self.totalBandpasses = []
        self.hardwareBandpasses = []

        componentList = ['detector.dat', 'm1.dat', 'm2.dat', 'm3.dat',
                         'lens1.dat', 'lens2.dat', 'lens3.dat']
        hardwareComponents = []
        for c in componentList:
            hardwareComponents.append(os.path.join(lsst.utils.getPackageDir('throughputs'),'baseline',c))

        self.bandpasses = ['u', 'g', 'r', 'i', 'z', 'y']
        for b in self.bandpasses:
            filterName = os.path.join(lsst.utils.getPackageDir('throughputs'),'baseline','filter_%s.dat' % b)
            components = hardwareComponents + [filterName]
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughputList(components)
            self.hardwareBandpasses.append(bandpassDummy)
            components = components + [os.path.join(lsst.utils.getPackageDir('throughputs'),'baseline','atmos.dat')]
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughputList(components)
            self.totalBandpasses.append(bandpassDummy)

    def testApplyIGM(self):

        """Test application of IGM from Lookup Tables to SED objects"""

        #Test that a warning comes up if input redshift is out of range and that no changes occurs to SED
        testSed = Sed()
        testSed.readSED_flambda(os.environ['SIMS_SED_LIBRARY_DIR'] + '/galaxySED/Inst.80E09.25Z.spec.gz')
        testFlambda = []
        for fVal in testSed.flambda:
            testFlambda.append(fVal)
        testIGM = ApplyIGM()
        testIGM.initializeIGM()
        with warnings.catch_warnings(record=True) as wa:
            testIGM.applyIGM(1.1, testSed)
            self.assertEqual(len(wa), 1)
            self.assertTrue('IGM Lookup tables' in str(wa[-1].message))
        np.testing.assert_equal(testFlambda, testSed.flambda)

        #Test that lookup table is read in correctly
        testTable15 = np.genfromtxt(str(os.environ['SIMS_SED_LIBRARY_DIR'] + '/igm/' +
                                        'MeanLookupTable_zSource1.5.tbl'))
        np.testing.assert_equal(testTable15, testIGM.meanLookups['1.5'])

        #Test output by making sure that an incoming sed with flambda = 1.0 everywhere will return the
        #transmission values of the lookup table as its flambda output
        testSed.setSED(testSed.wavelen, flambda = np.ones(len(testSed.wavelen)))
        testIGM.applyIGM(1.5, testSed)
        testTable15Above300 = testTable15[np.where(testTable15[:,0] >= 300.0)]
        testSed.resampleSED(wavelen_match = testTable15Above300[:,0])
        np.testing.assert_allclose(testTable15Above300[:,1], testSed.flambda, 1e-4)

    def loadmltSEDs(self, subset = None):

        """
        By default will load all seds in mlt directory. The user can also define a subset of
        what's in the directory and load only those SEDs instead. Will skip over extraneous
        files in sed folder.

        @param [in] subset is the list of the subset of files wanted if one doesn't want all files
        in the mlt directory.

        @param [out] sedList is the set of model SED spectra objects to be passed onto the matching
        routines.
        """

        files = []

        if subset is None:
            for fileName in os.listdir(self.mltDir):
                files.append(fileName)
        else:
            for fileName in subset:
                files.append(fileName)

        numFiles = len(files)
        numOn = 0

        sedList = []

        for fileName in files:
            if numOn % 100 == 0:
                print 'Loading %i of %i: MLT SEDs' % (numOn, numFiles)

            try:
                spec = Sed()
                spec.readSED_flambda(str(self.mltDir + '/' + fileName))
                spec.name = fileName

            except:
                continue

            sedList.append(spec)

            numOn += 1

        return sedList

    def testCalcBasicColors(self):

        """Tests the calculation of the colors of an SED in given bandpasses."""

        testUtils = matchBase()
        testSED = Sed()
        testPhot = BandpassDict.loadTotalBandpassesFromFiles(self.filterList,
                                        bandpassDir = os.path.join(lsst.utils.getPackageDir('throughputs'),'sdss'),
                                        bandpassRoot = 'sdss_')

        testSED.readSED_flambda(str(self.galDir + os.listdir(self.galDir)[0]))
        testMags = testPhot.magListForSed(testSED)
        testColors = []
        for filtNum in range(0, len(self.filterList)-1):
            testColors.append(testMags[filtNum] - testMags[filtNum+1])

        testOutput = testUtils.calcBasicColors([testSED], testPhot)
        np.testing.assert_equal([testColors], testOutput)

 def testPhotometricIndicesRaw(self):
     """
     Use manMagCalc_list with specified indices on an Sed.  Make sure
     that the appropriate magnitudes are or are not Nan
     """
     starName = os.path.join(getPackageDir('sims_sed_library'), defaultSpecMap['km20_5750.fits_g40_5790'])
     starPhot = BandpassDict.loadTotalBandpassesFromFiles()
     testSed = Sed()
     testSed.readSED_flambda(starName)
     indices = [1, 3]
     mags = starPhot.magListForSed(testSed, indices=indices)
     np.testing.assert_equal(mags[0], np.NaN)
     self.assertFalse(np.isnan(mags[1]), msg='mags[1] is NaN; should not be')
     np.testing.assert_equal(mags[2], np.NaN)
     self.assertFalse(np.isnan(mags[3]), msg='mags[3] is NaN; should not be')
     np.testing.assert_equal(mags[4], np.NaN)
     np.testing.assert_equal(mags[5], np.NaN)
     self.assertEqual(len(mags), 6)

    def testInitializeIGM(self):

        "Test Initialization Method"

        #Make sure that if we initialize IGM with new inputs that it is initializing with them
        testIGM = ApplyIGM()
        testSed = Sed()
        testSed.readSED_flambda(os.environ['SIMS_SED_LIBRARY_DIR'] + '/galaxySED/Inst.80E09.25Z.spec.gz')
        testIGM.applyIGM(1.8, testSed)
        testZmin = 1.8
        testZmax = 2.2
        #Want new values for testing,
        #so make sure we are not just putting in the same values as are already there
        self.assertNotEqual(testZmin, testIGM.zMin)
        self.assertNotEqual(testZmax, testIGM.zMax)
        testIGM.initializeIGM(zMin = testZmin, zMax = testZmax)
        self.assertEqual(testZmin, testIGM.zMin)
        self.assertEqual(testZmax, testIGM.zMax)