Python Sed.setSED方法代码示例

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    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

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    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)

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    def calcMagNorm(self, objectMags, sedObj, bandpassDict, mag_error = None,
                    redshift = None, filtRange = None):

        """
        This will find the magNorm value that gives the closest match to the magnitudes of the object
        using the matched SED. Uses scipy.optimize.leastsq to find the values of fluxNorm that minimizes
        the function: ((flux_obs - (fluxNorm*flux_model))/flux_error)**2.

        @param [in] objectMags are the magnitude values for the object with extinction matching that of
        the SED object. In the normal case using the selectSED routines above it will be dereddened mags.

        @param [in] sedObj is an Sed class instance that is set with the wavelength and flux of the
        matched SED

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

        @param [in] mag_error are provided error values for magnitudes in objectMags. If none provided
        then this defaults to 1.0. This should be an array of the same length as objectMags.

        @param [in] redshift is the redshift of the object if the magnitude is observed

        @param [in] filtRange is a selected range of filters specified by their indices in the bandpassList
        to match up against. Used when missing data in some magnitude bands.

        @param [out] bestMagNorm is the magnitude normalization for the given magnitudes and SED
        """

        import scipy.optimize as opt

        sedTest = Sed()
        sedTest.setSED(sedObj.wavelen, flambda = sedObj.flambda)
        if redshift is not None:
            sedTest.redshiftSED(redshift)
        imSimBand = Bandpass()
        imSimBand.imsimBandpass()
        zp = -2.5*np.log10(3631)  #Note using default AB zeropoint
        flux_obs = np.power(10,(objectMags + zp)/(-2.5))
        sedTest.resampleSED(wavelen_match=bandpassDict.values()[0].wavelen)
        sedTest.flambdaTofnu()
        flux_model = sedTest.manyFluxCalc(bandpassDict.phiArray, bandpassDict.wavelenStep)
        if filtRange is not None:
            flux_obs = flux_obs[filtRange]
            flux_model = flux_model[filtRange]
        if mag_error is None:
            flux_error = np.ones(len(flux_obs))
        else:
            flux_error = np.abs(flux_obs*(np.log(10)/(-2.5))*mag_error)
        bestFluxNorm = opt.leastsq(lambda x: ((flux_obs - (x*flux_model))/flux_error), 1.0)[0][0]
        sedTest.multiplyFluxNorm(bestFluxNorm)
        bestMagNorm = sedTest.calcMag(imSimBand)
        return bestMagNorm

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    def testMatchToRestFrame(self):
        """Test that Galaxies with no effects added into catalog mags are matched correctly."""
        np.random.seed(42)
        galPhot = BandpassDict.loadTotalBandpassesFromFiles()

        imSimBand = Bandpass()
        imSimBand.imsimBandpass()

        testMatching = selectGalaxySED(galDir = self.testSpecDir)
        testSEDList = testMatching.loadBC03()

        testSEDNames = []
        testMags = []
        testMagNormList = []
        magNormStep = 1

        for testSED in testSEDList:

            getSEDMags = Sed()
            testSEDNames.append(testSED.name)
            getSEDMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
            testMagNorm = np.round(np.random.uniform(20.0,22.0),magNormStep)
            testMagNormList.append(testMagNorm)
            fluxNorm = getSEDMags.calcFluxNorm(testMagNorm, imSimBand)
            getSEDMags.multiplyFluxNorm(fluxNorm)
            testMags.append(galPhot.magListForSed(getSEDMags))

        #Also testing to make sure passing in non-default bandpasses works
        #Substitute in nan values to simulate incomplete data.
        testMags[0][1] = np.nan
        testMags[0][2] = np.nan
        testMags[0][4] = np.nan
        testMags[1][1] = np.nan
        testMatchingResults = testMatching.matchToRestFrame(testSEDList, testMags,
                                                            bandpassDict = galPhot)
        self.assertEqual(None, testMatchingResults[0][0])
        self.assertEqual(testSEDNames[1:], testMatchingResults[0][1:])
        self.assertEqual(None, testMatchingResults[1][0])
        np.testing.assert_almost_equal(testMagNormList[1:], testMatchingResults[1][1:], decimal = magNormStep)

        #Test Match Errors
        errMags = np.array((testMags[2], testMags[2], testMags[2], testMags[2]))
        errMags[1,1] += 1. #Total MSE will be 2/(5 colors) = 0.4
        errMags[2, 0:2] = np.nan
        errMags[2, 3] += 1. #Total MSE will be 2/(3 colors) = 0.667
        errMags[3, :] = None
        errSED = testSEDList[2]
        testMatchingResultsErrors = testMatching.matchToRestFrame([errSED], errMags,
                                                                  bandpassDict = galPhot)
        np.testing.assert_almost_equal(np.array((0.0, 0.4, 2./3.)), testMatchingResultsErrors[2][0:3],
                                       decimal = 3)
        self.assertEqual(None, testMatchingResultsErrors[2][3])

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    def testSEDCopyBasicColors(self):

        """Tests that when makeCopy=True in calcBasicColors the SED object is unchanged after calling
        and that colors are still accurately calculated"""

        testUtils = matchBase()
        testSED = Sed()
        copyTest = 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]))
        copyTest.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
        testLambda = copyTest.wavelen[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([copyTest], testPhot, makeCopy=True)

        self.assertEqual(testLambda, copyTest.wavelen[0])
        np.testing.assert_equal([testColors], testOutput)

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    def testCalcMagNorm(self):

        """Tests the calculation of magnitude normalization for an SED with the given magnitudes
        in the given bandpasses."""

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

        unChangedSED = Sed()
        unChangedSED.readSED_flambda(str(self.galDir + os.listdir(self.galDir)[0]))

        imSimBand = Bandpass()
        imSimBand.imsimBandpass()
        testSED = Sed()
        testSED.setSED(unChangedSED.wavelen, flambda = unChangedSED.flambda)
        magNorm = 20.0
        redVal = 0.1
        testSED.redshiftSED(redVal)
        fluxNorm = testSED.calcFluxNorm(magNorm, imSimBand)
        testSED.multiplyFluxNorm(fluxNorm)
        sedMags = testPhot.magListForSed(testSED)
        stepSize = 0.001
        testMagNorm = testUtils.calcMagNorm(sedMags, unChangedSED, testPhot, redshift = redVal)
        # Test adding in mag_errors. If an array of np.ones is passed in we should get same result
        testMagNormWithErr = testUtils.calcMagNorm(sedMags, unChangedSED, testPhot,
                                                   mag_error = np.ones(len(sedMags)), redshift = redVal)
        # Also need to add in test for filtRange
        sedMagsIncomp = sedMags
        sedMagsIncomp[1] = None
        filtRangeTest = [0, 2, 3, 4]
        testMagNormFiltRange = testUtils.calcMagNorm(sedMagsIncomp, unChangedSED, testPhot,
                                                     redshift = redVal, filtRange = filtRangeTest)
        self.assertAlmostEqual(magNorm, testMagNorm, delta = stepSize)
        self.assertAlmostEqual(magNorm, testMagNormWithErr, delta = stepSize)
        self.assertAlmostEqual(magNorm, testMagNormFiltRange, delta = stepSize)

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    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()
        sedName = os.path.join(getPackageDir('sims_sed_library'), 'galaxySED')
        testSed.readSED_flambda(os.path.join(sedName,
                                             'Burst.10E08.002Z.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.assertIn('IGM Lookup tables', str(wa[-1].message))
        np.testing.assert_equal(testFlambda, testSed.flambda)

        # Test that lookup table is read in correctly
        testTable15 = np.genfromtxt(str(getPackageDir('sims_catUtils') +
                                        '/python/lsst/sims/catUtils/IGM/igm_tables/' +
                                        'MeanLookupTable_zSource1.5.tbl.gz'))
        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)

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]

#.........这里部分代码省略.........
                """
                Info about the following population cuts:
                From Zeljko: "This color corresponds to the temperature (roughly spectral type M0) where
                Kurucz models become increasingly bad, and thus we switch to empirical SEDs (the problem
                is that for M and later stars, the effective surface temperature is low enough for
                molecules to form, and their opacity is too complex to easily model, especially TiO)."
                """
                mIn = np.where(((pop < 10) | (pop >= 20)) & (sDSSunred[:,2] - sDSSunred[:,3] > 0.59))
                kIn = np.where(((pop < 10) | (pop >= 20)) & (sDSSunred[:,2] - sDSSunred[:,3] <= 0.59))
                hIn = np.where((pop >= 10) & (pop < 15))
                heIn = np.where((pop >= 15) & (pop < 20))

                sEDNameK, magNormK, matchErrorK = selectStarSED0.findSED(listDict['kurucz'],
                                                                         sDSSunred[kIn], ra[kIn], dec[kIn],
                                                                         reddening = False,
                                                                         colors = colorDict['kurucz'])
                sEDNameM, magNormM, matchErrorM = selectStarSED0.findSED(listDict['mlt'],
                                                                         sDSSunred[mIn], ra[mIn], dec[mIn],
                                                                         reddening = False,
                                                                         colors = colorDict['mlt'])
                sEDNameH, magNormH, matchErrorH = selectStarSED0.findSED(listDict['H'],
                                                                         sDSSunred[hIn], ra[hIn], dec[hIn],
                                                                         reddening = False,
                                                                         colors = colorDict['H'])
                sEDNameHE, magNormHE, matchErrorHE = selectStarSED0.findSED(listDict['HE'],
                                                                            sDSSunred[heIn],
                                                                            ra[heIn], dec[heIn],
                                                                            reddening = False,
                                                                            colors = colorDict['HE'])
                chunkNames = np.empty(readSize, dtype = 'S32')
                chunkTypes = np.empty(readSize, dtype = 'S8')
                chunkMagNorms = np.zeros(readSize)
                chunkMatchErrors = np.zeros(readSize)
                chunkNames[kIn] = sEDNameK
                chunkTypes[kIn] = 'kurucz'
                chunkMagNorms[kIn] = magNormK
                chunkMatchErrors[kIn] = matchErrorK
                chunkNames[mIn] = sEDNameM
                chunkTypes[mIn] = 'mlt'
                chunkMagNorms[mIn] = magNormM
                chunkMatchErrors[mIn] = matchErrorM
                chunkNames[hIn] = sEDNameH
                chunkTypes[hIn] = 'H'
                chunkMagNorms[hIn] = magNormH
                chunkMatchErrors[hIn] = matchErrorH
                chunkNames[heIn] = sEDNameHE
                chunkTypes[heIn] = 'HE'
                chunkMagNorms[heIn] = magNormHE
                chunkMatchErrors[heIn] = matchErrorHE
                lsstMagsUnred = []
                for sedName, sedType, magNorm, matchError in zip(chunkNames.astype(str),
                                                                 chunkTypes.astype(str),
                                                                 chunkMagNorms,
                                                                 chunkMatchErrors):
                    testSED = Sed()
                    testSED.setSED(listDict[sedType][positionDict[sedName]].wavelen,
                                   flambda = listDict[sedType][positionDict[sedName]].flambda)
                    fluxNorm = testSED.calcFluxNorm(magNorm, imSimBand)
                    testSED.multiplyFluxNorm(fluxNorm)
                    lsstMagsUnred.append(lsstPhot.magListForSed(testSED))
                #If the extinction value is negative then it will add the reddening back in
                lsstMags = selectStarSED0.deReddenMags((-1.0*am), lsstMagsUnred,
                                                       lsstExtCoeffs)
                distKpc = self.convDMtoKpc(DM)
                ebv = am / 2.285 #From Schlafly and Finkbeiner 2011, (ApJ, 737, 103)  for sdssr
                ebvInf = amInf / 2.285
                for line in range(0, readSize):
                    outFmt = '%i,%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,' +\
                             '%3.7f,%3.7f,%s,' +\
                             '%3.7f,%3.7f,' +\
                             '%3.7f,%3.7f,%3.7f,' +\
                             '%3.7f,%3.7f,%3.7f,' +\
                             '%3.7f,%3.7f,%3.7f,%3.7f,' +\
                             '%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,' +\
                             '%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,' +\
                             '%3.7f,%i,%3.7f,%3.7f,%3.7f\n'
                    if readSize == 1:
                        if inFits == True:
                            sDSS = sDSS[0]
                        outDat = (oID, ra[line], dec[line], gall, galb, coordX,
                                  coordY, coordZ, chunkNames,
                                  chunkMagNorms, chunkMatchErrors,
                                  lsstMags[line][0], lsstMags[line][1], lsstMags[line][2],
                                  lsstMags[line][3], lsstMags[line][4], lsstMags[line][5],
                                  sDSS[0], sDSS[1], sDSS[2], sDSS[3],
                                  sDSS[4], absSDSSr, pmRA, pmDec, vRad,
                                  pml, pmb, vRadlb, vR, vPhi, vZ,
                                  FeH, pop, distKpc, ebv, ebvInf)
                    else:
                        outDat = (oID[line], ra[line], dec[line], gall[line], galb[line], coordX[line],
                                  coordY[line], coordZ[line], chunkNames[line],
                                  chunkMagNorms[line], chunkMatchErrors[line],
                                  lsstMags[line][0], lsstMags[line][1], lsstMags[line][2],
                                  lsstMags[line][3], lsstMags[line][4], lsstMags[line][5],
                                  sDSS[line][0], sDSS[line][1], sDSS[line][2], sDSS[line][3],
                                  sDSS[line][4], absSDSSr[line], pmRA[line], pmDec[line], vRad[line],
                                  pml[line], pmb[line], vRadlb[line], vR[line], vPhi[line], vZ[line],
                                  FeH[line], pop[line], distKpc[line], ebv[line], ebvInf[line])
                    fOut.write(outFmt % outDat)
                print('Chunk Num Done = %i out of %i' % (chunk+1, numChunks))

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]

#.........这里部分代码省略.........
        @param [out] magNormMatches are the magnitude normalizations for the given magnitudes and
        matched SED.

        @param [out] matchErrors contains the Mean Squared Error between the colors of each object and 
        the colors of the matched SED.
        """

        #Set up photometry to calculate model Mags
        if bandpassDict is None:
            galPhot = BandpassDict.loadTotalBandpassesFromFiles(['u','g','r','i','z'],
                                            bandpassDir = os.path.join(lsst.utils.getPackageDir('throughputs'),'sdss'),
                                            bandpassRoot = 'sdss_')
        else:
            galPhot = bandpassDict

        #Calculate ebv from ra, dec coordinates if needed
        if reddening == True:
            #Check that catRA and catDec are included
            if catRA is None or catDec is None:
                raise RuntimeError("Reddening is True, but catRA and catDec are not included.")
            calcEBV = ebv()
            raDec = np.array((catRA,catDec))
            #If only matching one object need to reshape for calculateEbv
            if len(raDec.shape) == 1:
                raDec = raDec.reshape((2,1))
            ebvVals = calcEBV.calculateEbv(equatorialCoordinates = raDec)
            objMags = self.deReddenMags(ebvVals, catMags, extCoeffs)
        else:
            objMags = catMags

        minRedshift = np.round(np.min(catRedshifts), dzAcc)
        maxRedshift = np.round(np.max(catRedshifts), dzAcc)
        dz = np.power(10., (-1*dzAcc))

        redshiftRange = np.round(np.arange(minRedshift - dz, maxRedshift + (2*dz), dz), dzAcc)
        numRedshifted = 0
        sedMatches = [None] * len(catRedshifts)
        magNormMatches = [None] * len(catRedshifts)
        matchErrors = [None] * len(catRedshifts)
        redshiftIndex = np.argsort(catRedshifts)

        numOn = 0
        notMatched = 0
        lastRedshift = -100
        print 'Starting Matching. Arranged by redshift value.'
        for redshift in redshiftRange:

            if numRedshifted % 10 == 0:
                print '%i out of %i redshifts gone through' % (numRedshifted, len(redshiftRange))
            numRedshifted += 1

            colorSet = []
            for galSpec in sedList:
                sedColors = []
                fileSED = Sed()
                fileSED.setSED(wavelen = galSpec.wavelen, flambda = galSpec.flambda)
                fileSED.redshiftSED(redshift)
                sedColors = self.calcBasicColors([fileSED], galPhot, makeCopy = True)
                colorSet.append(sedColors)
            colorSet = np.transpose(colorSet)
            for currentIndex in redshiftIndex[numOn:]:
                matchMags = objMags[currentIndex]
                if lastRedshift < np.round(catRedshifts[currentIndex],dzAcc) <= redshift:
                    colorRange = np.arange(0, len(galPhot)-1)
                    matchColors = []
                    for colorNum in colorRange:
                        matchColors.append(matchMags[colorNum] - matchMags[colorNum+1])
                    #This is done to handle objects with incomplete magnitude data
                    filtNums = np.arange(0, len(galPhot))
                    if np.isnan(np.amin(matchColors))==True:
                        colorRange = np.where(np.isnan(matchColors)==False)[0]
                        filtNums = np.unique([colorRange, colorRange+1]) #Pick right filters in calcMagNorm
                    if len(colorRange) == 0:
                        print 'Could not match object #%i. No magnitudes for two adjacent bandpasses.' \
                              % (currentIndex)
                        notMatched += 1
                        #Don't need to assign 'None' here in result array, b/c 'None' is default value
                    else:
                        distanceArray = [np.zeros(len(sedList))]
                        for colorNum in colorRange:
                            distanceArray += np.power((colorSet[colorNum] - matchColors[colorNum]),2)
                        matchedSEDNum = np.nanargmin(distanceArray)
                        sedMatches[currentIndex] = sedList[matchedSEDNum].name
                        magNormVal = self.calcMagNorm(np.array(matchMags), sedList[matchedSEDNum], 
                                                      galPhot, mag_error = mag_error,
                                                      redshift = catRedshifts[currentIndex],
                                                      filtRange = filtNums)
                        magNormMatches[currentIndex] = magNormVal
                        matchErrors[currentIndex] = (distanceArray[0,matchedSEDNum]/len(colorRange))
                    numOn += 1
                else:
                    break
            lastRedshift = redshift

        print 'Done Matching. Matched %i of %i catalog objects to SEDs' % (len(catMags)-notMatched, 
                                                                           len(catMags))
        if notMatched > 0:
            print '%i objects did not get matched.' % (notMatched)

        return sedMatches, magNormMatches, matchErrors

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    def testFindSED(self):
        """Pull SEDs from each type and make sure that each SED gets matched to itself.
        Includes testing with extinction and passing in only colors."""
        np.random.seed(42)
        starPhot = BandpassDict.loadTotalBandpassesFromFiles(('u','g','r','i','z'),
                                        bandpassDir = os.path.join(lsst.utils.getPackageDir('throughputs'),'sdss'),
                                        bandpassRoot = 'sdss_')

        imSimBand = Bandpass()
        imSimBand.imsimBandpass()

        testMatching = selectStarSED(sEDDir = self.testSpecDir, kuruczDir = self.testKDir,
                                     mltDir = self.testMLTDir, wdDir = self.testWDDir)
        testSEDList = []
        testSEDList.append(testMatching.loadKuruczSEDs())
        testSEDList.append(testMatching.loadmltSEDs())
        testSEDListH, testSEDListHE = testMatching.loadwdSEDs()
        testSEDList.append(testSEDListH)
        testSEDList.append(testSEDListHE)

        testSEDNames = []
        testMags = []
        testMagNormList = []
        magNormStep = 1

        for typeList in testSEDList:
            if len(typeList) != 0:
                typeSEDNames = []
                typeMags = []
                typeMagNorms = []
                for testSED in typeList:
                    getSEDMags = Sed()
                    typeSEDNames.append(testSED.name)
                    getSEDMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
                    testMagNorm = np.round(np.random.uniform(20.0,22.0),magNormStep)
                    typeMagNorms.append(testMagNorm)
                    fluxNorm = getSEDMags.calcFluxNorm(testMagNorm, imSimBand)
                    getSEDMags.multiplyFluxNorm(fluxNorm)
                    typeMags.append(starPhot.magListForSed(getSEDMags))
                testSEDNames.append(typeSEDNames)
                testMags.append(typeMags)
                testMagNormList.append(typeMagNorms)

        fakeRA = np.ones(len(testSEDList[0]))
        fakeDec = np.ones(len(testSEDList[0]))

        #Since default bandpassDict should be SDSS ugrizy shouldn't need to specify it
        #Substitute in nan values to simulate incomplete data.
        for typeList, names, mags, magNorms in zip(testSEDList, testSEDNames, testMags, testMagNormList):
            if len(typeList) > 2:
                nanMags = np.array(mags)
                nanMags[0][0] = np.nan
                nanMags[0][2] = np.nan
                nanMags[0][3] = np.nan
                nanMags[1][1] = np.nan
                testMatchingResults = testMatching.findSED(typeList, nanMags, reddening = False)
                self.assertEqual(None, testMatchingResults[0][0])
                self.assertEqual(names[1:], testMatchingResults[0][1:])
                self.assertEqual(None, testMatchingResults[1][0])
                np.testing.assert_almost_equal(magNorms[1:], testMatchingResults[1][1:],
                                               decimal = magNormStep)
            else:
                testMatchingResults = testMatching.findSED(typeList, mags, reddening = False)
                self.assertEqual(names, testMatchingResults[0])
                np.testing.assert_almost_equal(magNorms, testMatchingResults[1], decimal = magNormStep)

        #Test Null Values option
        nullMags = np.array(testMags[0])
        nullMags[0][0] = -99.
        nullMags[0][4] = -99.
        nullMags[1][0] = -99.
        nullMags[1][1] = -99.
        testMatchingResultsNull = testMatching.findSED(testSEDList[0], nullMags,
                                                       nullValues = -99., reddening = False)
        self.assertEqual(testSEDNames[0], testMatchingResultsNull[0])
        np.testing.assert_almost_equal(testMagNormList[0], testMatchingResultsNull[1],
                                       decimal = magNormStep)

        #Test Error Output
        errMags = np.array((testMags[0][0], testMags[0][0], testMags[0][0], testMags[0][0]))
        errMags[1,1] += 1. #Total MSE will be 2/(4 colors) = 0.5
        errMags[2, 0:2] = np.nan
        errMags[2, 3] += 1. #Total MSE will be 2/(2 colors) = 1.0
        errMags[3, :] = None
        errSED = testSEDList[0][0]
        testMatchingResultsErrors = testMatching.findSED([errSED], errMags, reddening = False)
        np.testing.assert_almost_equal(np.array((0.0, 0.5, 1.0)), testMatchingResultsErrors[2][0:3],
                                       decimal = 3)
        self.assertEqual(None, testMatchingResultsErrors[2][3])

        #Now test what happens if we pass in a bandpassDict
        testMatchingResultsNoDefault = testMatching.findSED(testSEDList[0], testMags[0],
                                                            bandpassDict = starPhot,
                                                            reddening = False)
        self.assertEqual(testSEDNames[0], testMatchingResultsNoDefault[0])
        np.testing.assert_almost_equal(testMagNormList[0], testMatchingResultsNoDefault[1],
                                       decimal = magNormStep)

        #Test Reddening
        testRA = np.random.uniform(10,170,len(testSEDList[0]))
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    def testMatchToObserved(self):
        """Test that Galaxy SEDs with extinction or redshift are matched correctly"""
        np.random.seed(42)
        galPhot = BandpassDict.loadTotalBandpassesFromFiles()

        imSimBand = Bandpass()
        imSimBand.imsimBandpass()

        testMatching = selectGalaxySED(galDir = self.testSpecDir)
        testSEDList = testMatching.loadBC03()

        testSEDNames = []
        testRA = []
        testDec = []
        testRedshifts = []
        testMagNormList = []
        magNormStep = 1
        extCoeffs = [1.8140, 1.4166, 0.9947, 0.7370, 0.5790, 0.4761]
        testMags = []
        testMagsRedshift = []
        testMagsExt = []

        for testSED in testSEDList:

            #As a check make sure that it matches when no extinction and no redshift are present
            getSEDMags = Sed()
            testSEDNames.append(testSED.name)
            getSEDMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
            testMags.append(galPhot.magListForSed(getSEDMags))

            #Check Extinction corrections
            sedRA = np.random.uniform(10,170)
            sedDec = np.random.uniform(10,80)
            testRA.append(sedRA)
            testDec.append(sedDec)
            raDec = np.array((sedRA, sedDec)).reshape((2,1))
            ebvVal = ebv().calculateEbv(equatorialCoordinates = raDec)
            extVal = ebvVal*extCoeffs
            testMagsExt.append(galPhot.magListForSed(getSEDMags) + extVal)

            #Setup magnitudes for testing matching to redshifted values
            getRedshiftMags = Sed()
            testZ = np.round(np.random.uniform(1.1,1.3),3)
            testRedshifts.append(testZ)
            testMagNorm = np.round(np.random.uniform(20.0,22.0),magNormStep)
            testMagNormList.append(testMagNorm)
            getRedshiftMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
            getRedshiftMags.redshiftSED(testZ)
            fluxNorm = getRedshiftMags.calcFluxNorm(testMagNorm, imSimBand)
            getRedshiftMags.multiplyFluxNorm(fluxNorm)
            testMagsRedshift.append(galPhot.magListForSed(getRedshiftMags))

        #Will also test in passing of non-default bandpass
        testNoExtNoRedshift = testMatching.matchToObserved(testSEDList, testMags, np.zeros(20),
                                                           reddening = False,
                                                           bandpassDict = galPhot)
        testMatchingEbvVals = testMatching.matchToObserved(testSEDList, testMagsExt, np.zeros(20),
                                                           catRA = testRA, catDec = testDec,
                                                           reddening = True, extCoeffs = extCoeffs,
                                                           bandpassDict = galPhot)
        #Substitute in nan values to simulate incomplete data and make sure magnorm works too.
        testMagsRedshift[0][1] = np.nan
        testMagsRedshift[0][3] = np.nan
        testMagsRedshift[0][4] = np.nan
        testMagsRedshift[1][1] = np.nan
        testMatchingRedshift = testMatching.matchToObserved(testSEDList, testMagsRedshift, testRedshifts,
                                                            dzAcc = 3, reddening = False,
                                                            bandpassDict = galPhot)

        self.assertEqual(testSEDNames, testNoExtNoRedshift[0])
        self.assertEqual(testSEDNames, testMatchingEbvVals[0])
        self.assertEqual(None, testMatchingRedshift[0][0])
        self.assertEqual(testSEDNames[1:], testMatchingRedshift[0][1:])
        self.assertEqual(None, testMatchingRedshift[1][0])
        np.testing.assert_almost_equal(testMagNormList[1:], testMatchingRedshift[1][1:],
                                       decimal = magNormStep)

        #Test Match Errors
        errMag = testMagsRedshift[2]
        errRedshift = testRedshifts[2]
        errMags = np.array((errMag, errMag, errMag, errMag))
        errRedshifts = np.array((errRedshift, errRedshift, errRedshift, errRedshift))
        errMags[1,1] += 1. #Total MSE will be 2/(5 colors) = 0.4
        errMags[2, 0:2] = np.nan
        errMags[2, 3] += 1. #Total MSE will be 2/(3 colors) = 0.667
        errMags[3, :] = None
        errSED = testSEDList[2]
        testMatchingResultsErrors = testMatching.matchToObserved([errSED], errMags, errRedshifts,
                                                                 reddening = False,
                                                                 bandpassDict = galPhot,
                                                                 dzAcc = 3)
        np.testing.assert_almost_equal(np.array((0.0, 0.4, 2./3.)), testMatchingResultsErrors[2][0:3],
                                       decimal = 2) #Give a little more leeway due to redshifting effects
        self.assertEqual(None, testMatchingResultsErrors[2][3])

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setSED [as 别名]
    def applyIGM(self, redshift, sedobj):

        """
        Apply IGM extinction to already redshifted sed with redshift  
        between zMin and zMax defined by range of lookup tables
        
        @param [in] redshift is the redshift of the incoming SED object

        @param [in] sedobj is the SED object to which IGM extinction will be applied. This object
        will be modified as a result of this.
        """

        if self.IGMisInitialized == False:
            self.initializeIGM()

        # First make sure redshift is in range of lookup tables.
        if (redshift < self.zMin) or (redshift > self.zMax):
            warnings.warn(
                str(
                    "IGM Lookup tables only applicable for "
                    + str(self.zMin)
                    + " < z < "
                    + str(self.zMax)
                    + ". No action taken"
                )
            )
            return

        # Now read in closest two lookup tables for given redshift
        lowerSed = Sed()
        upperSed = Sed()
        for lower, upper in zip(self.zRange[:-1], self.zRange[1:]):
            if lower <= redshift <= upper:
                lowerSed.setSED(self.meanLookups[str(lower)][:, 0], flambda=self.meanLookups[str(lower)][:, 1])
                upperSed.setSED(self.meanLookups[str(upper)][:, 0], flambda=self.meanLookups[str(upper)][:, 1])
                break

        # Redshift lookup tables to redshift of source, i.e. if source redshift is 1.78 shift lookup
        # table for 1.7 and lookup table for 1.8 to up and down to 1.78, respectively
        zLowerShift = ((1.0 + redshift) / (1.0 + lower)) - 1.0
        zUpperShift = ((1.0 + redshift) / (1.0 + upper)) - 1.0
        lowerSed.redshiftSED(zLowerShift)
        upperSed.redshiftSED(zUpperShift)

        # Resample lower and upper transmission data onto same wavelength grid.
        minWavelen = 300.0  # All lookup tables are usable above 300nm
        maxWavelen = np.amin([lowerSed.wavelen[-1], upperSed.wavelen[-1]]) - 0.01
        lowerSed.resampleSED(wavelen_min=minWavelen, wavelen_max=maxWavelen, wavelen_step=0.01)
        upperSed.resampleSED(wavelen_match=lowerSed.wavelen)

        # Now insert this into a transmission array of 1.0 beyond the limits of current application
        # So that we can get an sed back that extends to the longest wavelengths of the incoming sed
        finalWavelen = np.arange(300.0, sedobj.wavelen[-1] + 0.01, 0.01)
        finalFlambdaExtended = np.ones(len(finalWavelen))

        # Weighted Average of Transmission from each lookup table to get final transmission
        # table at desired redshift
        dzGrid = self.zDelta  # Step in redshift between transmission lookup table files
        finalSed = Sed()
        finalFlambda = lowerSed.flambda * (1.0 - ((redshift - lower) / dzGrid)) + upperSed.flambda * (
            1.0 - ((upper - redshift) / dzGrid)
        )
        finalFlambdaExtended[0 : len(finalFlambda)] = finalFlambda
        finalSed.setSED(wavelen=finalWavelen, flambda=finalFlambdaExtended)

        # Resample incoming sed to new grid so that we don't get warnings from multiplySED
        # about matching wavelength grids
        sedobj.resampleSED(wavelen_match=finalSed.wavelen)

        # Now multiply transmission curve by input SED to get final result and make it the new flambda
        # data in the original sed which also is now on a new grid starting at 300 nm
        test = sedobj.multiplySED(finalSed)
        sedobj.flambda = test.flambda