Python Sed.calcFluxNorm方法代码示例

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

        bandpassDir = os.path.join(lsst.utils.getPackageDir("sims_photUtils"), "tests", "cartoonSedTestData")

        cartoon_dict = BandpassDict.loadTotalBandpassesFromFiles(
            ["u", "g", "r", "i", "z"], bandpassDir=bandpassDir, bandpassRoot="test_bandpass_"
        )

        testBandPasses = {}
        keys = ["u", "g", "r", "i", "z"]

        bplist = []

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

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

        sedFileName = os.path.join(lsst.utils.getPackageDir("sims_sed_library"), "starSED", "kurucz")
        sedFileName = os.path.join(sedFileName, "km20_5750.fits_g40_5790.gz")
        ss = Sed()
        ss.readSED_flambda(sedFileName)

        controlBandpass = Bandpass()
        controlBandpass.imsimBandpass()
        ff = ss.calcFluxNorm(22.0, controlBandpass)
        ss.multiplyFluxNorm(ff)

        testMags = cartoon_dict.magListForSed(ss)

        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(testMags))
        self.assertTrue(len(mags) > 0)
        for j in range(len(mags)):
            self.assertAlmostEqual(mags[j], testMags[j], 10)

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

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import calcFluxNorm [as 别名]
class uncertaintyUnitTest(unittest.TestCase):
    """
    Test the calculation of photometric uncertainties
    """

    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 tearDown(self):
        del self.starSED
        del self.bandpasses
        del self.hardwareBandpasses
        del self.totalBandpasses

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import calcFluxNorm [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 calcFluxNorm [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 calcFluxNorm [as 别名]
    def testAlternateBandpassesGalaxies(self):
        """
        the same as testAlternateBandpassesStars, but for galaxies
        """

        obs_metadata_pointed = ObservationMetaData(mjd=50000.0,
                                                   boundType='circle',
                                                   pointingRA=0.0, pointingDec=0.0,
                                                   boundLength=10.0)

        dtype = np.dtype([('galid', np.int),
                          ('ra', np.float),
                          ('dec', np.float),
                          ('uTotal', np.float),
                          ('gTotal', np.float),
                          ('rTotal', np.float),
                          ('iTotal', np.float),
                          ('zTotal', np.float),
                          ('uBulge', np.float),
                          ('gBulge', np.float),
                          ('rBulge', np.float),
                          ('iBulge', np.float),
                          ('zBulge', np.float),
                          ('uDisk', np.float),
                          ('gDisk', np.float),
                          ('rDisk', np.float),
                          ('iDisk', np.float),
                          ('zDisk', np.float),
                          ('uAgn', np.float),
                          ('gAgn', np.float),
                          ('rAgn', np.float),
                          ('iAgn', np.float),
                          ('zAgn', np.float),
                          ('bulgeName', str, 200),
                          ('bulgeNorm', np.float),
                          ('bulgeAv', np.float),
                          ('diskName', str, 200),
                          ('diskNorm', np.float),
                          ('diskAv', np.float),
                          ('agnName', str, 200),
                          ('agnNorm', np.float),
                          ('redshift', np.float)])

        test_cat = cartoonGalaxies(self.galaxy, obs_metadata=obs_metadata_pointed)
        with lsst.utils.tests.getTempFilePath('.txt') as catName:
            test_cat.write_catalog(catName)
            catData = np.genfromtxt(catName, dtype=dtype, delimiter=', ')

        self.assertGreater(len(catData), 0)

        cartoonDir = getPackageDir('sims_photUtils')
        cartoonDir = os.path.join(cartoonDir, 'tests', 'cartoonSedTestData')
        sedDir = getPackageDir('sims_sed_library')

        testBandpasses = {}
        keys = ['u', 'g', 'r', 'i', 'z']

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

        imsimBand = Bandpass()
        imsimBand.imsimBandpass()

        specMap = defaultSpecMap

        ct = 0
        for line in catData:
            bulgeMagList = []
            diskMagList = []
            agnMagList = []
            if line['bulgeName'] == 'None':
                for bp in keys:
                    np.testing.assert_equal(line['%sBulge' % bp], np.NaN)
                    bulgeMagList.append(np.NaN)
            else:
                ct += 1
                dummySed = Sed()
                dummySed.readSED_flambda(os.path.join(sedDir, specMap[line['bulgeName']]))
                fnorm = dummySed.calcFluxNorm(line['bulgeNorm'], imsimBand)
                dummySed.multiplyFluxNorm(fnorm)
                a_int, b_int = dummySed.setupCCM_ab()
                dummySed.addDust(a_int, b_int, A_v=line['bulgeAv'])
                dummySed.redshiftSED(line['redshift'], dimming=True)
                dummySed.resampleSED(wavelen_match=testBandpasses['u'].wavelen)
                for bpName in keys:
                    mag = dummySed.calcMag(testBandpasses[bpName])
                    self.assertAlmostEqual(mag, line['%sBulge' % bpName], 10)
                    bulgeMagList.append(mag)

            if line['diskName'] == 'None':
                for bp in keys:
                    np.assert_equal(line['%sDisk' % bp], np.NaN)
                    diskMagList.append(np.NaN)
            else:
                ct += 1
                dummySed = Sed()
                dummySed.readSED_flambda(os.path.join(sedDir, specMap[line['diskName']]))
                fnorm = dummySed.calcFluxNorm(line['diskNorm'], imsimBand)
                dummySed.multiplyFluxNorm(fnorm)
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import calcFluxNorm [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 calcFluxNorm [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 calcFluxNorm [as 别名]
class uncertaintyUnitTest(unittest.TestCase):
    """
    Test the calculation of photometric uncertainties
    """

    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 tearDown(self):
        del self.starSED
        del self.bandpasses
        del self.hardwareBandpasses
        del self.totalBandpasses

    def testUncertaintyExceptions(self):
        """
        Test that calcSNR_m5 raises exceptions when it needs to
        """
        totalDict, hardwareDict = BandpassDict.loadBandpassesFromFiles()
        magnitudes = numpy.array([22.0, 23.0, 24.0, 25.0, 26.0, 27.0])
        shortMagnitudes = numpy.array([22.0])
        photParams = PhotometricParameters()
        shortGamma = numpy.array([1.0, 1.0])
        self.assertRaises(RuntimeError, calcSNR_m5, magnitudes, totalDict.values(), shortMagnitudes, photParams)
        self.assertRaises(RuntimeError, calcSNR_m5, shortMagnitudes, totalDict.values(), magnitudes, photParams)
        self.assertRaises(
            RuntimeError, calcSNR_m5, magnitudes, totalDict.values(), magnitudes, photParams, gamma=shortGamma
        )
        snr, gg = calcSNR_m5(magnitudes, totalDict.values(), magnitudes, photParams)

    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
            )
#.........这里部分代码省略.........