Python Planck13.distmod方法代碼示例

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import distmod [as 別名]
    def mkParamDict(self,zcontrol):
        """Make a dictionary to store all info about parameters"""
        from txtobj import txtobj
        pf = txtobj(self.options.mcmcparamfile)

        if self.options.twogauss:
            pf.use[pf.param == 'popB2mean'] = 1
            pf.use[pf.param == 'popB2std'] = 1
        if self.options.skewedgauss:
            pf.use[pf.param == 'skewB'] = 1
        if len(self.options.fix):
            for fixvar in self.options.fix:
                print('Fixing parameter %s!!'%fixvar)
                pf.fixed[pf.param == fixvar] = 1
        if len(self.options.use):
            for use in self.options.use:
                usevar,useval = use
                useval = int(useval)
                print('use = %i for parameter %s!!'%(useval,usevar))
                pf.use[pf.param == usevar] = useval
        if len(self.options.bounds):
            for bounds in self.options.bounds:
                boundsvar,lbound,ubound = bounds
                lbound,ubound = float(lbound),float(ubound)
                print('%.3f < %s < %.3f !!'%(lbound,boundsvar,ubound))
                pf.lbound[pf.param == boundsvar] = lbound
                pf.ubound[pf.param == boundsvar] = ubound
        if len(self.options.guess):
            for guess in self.options.guess:
                guessvar,guessval = guess
                guessvar = float(guessvar)
                print('initial guess = %.3f for parameter %s!!'%(
                        guessval,guessvar))
                pf.guess[pf.param == guessvar] = guessval
        if len(self.options.prior):
            for prior in self.options.prior:
                priorvar,priormean,priorstd = prior
                priormean,priorstd = float(priormean),float(priorstd)
                print('Prior = %.3f +/- %.3f for parameter %s!!'%(
                        priormean,priorstd,priorvar))
                pf.prior[pf.param == priorvar] = priormean
                pf.sigma[pf.param == priorvar] = sigma
        if len(self.options.bins):
            for bins in self.options.bins:
                binvar,nbins = bins
                nbins = float(nbins)
                print('%i bins for parameter %s!!'%(nbins,binvar))
                pf.bins[pf.param == binvar] = nbins


        self.pardict = {}
        idx = 0
        for par,i in zip(pf.param,xrange(len(pf.param))):
            self.pardict[par] = {'guess':pf.guess[i],'prior_mean':pf.prior[i],
                                 'prior_std':pf.sigma[i],'fixed':pf.fixed[i],
                                 'use':pf.use[i],'addcosmo':pf.addcosmo[i],
                                 'mcstep':self.options.mcrandstep,
                                 'bins':pf.bins[i],'zpoly':pf.zpoly[i],
                                 'bounds':(pf.lbound[i],pf.ubound[i])}
            if not pf.use[i]: self.pardict[par]['idx'] = -1
            else:
                self.pardict[par]['idx'] = idx
                if pf.addcosmo[i]:
                    if pf.bins[i] == 1:
                        self.pardict[par]['guess'] = pf.guess[i] + cosmo.distmod(zcontrol).value
                        self.pardict[par]['prior_mean'] = pf.prior[i] + cosmo.distmod(zcontrol).value
                        self.pardict[par]['idx'] = idx + np.arange(len(zcontrol))
                        self.pardict[par]['bounds'] = (pf.lbound[i] + cosmo.distmod(zcontrol).value,
                                                       pf.ubound[i] + cosmo.distmod(zcontrol).value)
                        if pf.use[i]: idx += len(zcontrol)
                    else:
                        zcontrolCC = np.logspace(np.log10(min(zcontrol)),np.log10(max(zcontrol)),pf.bins[i])
                        self.pardict[par]['guess'] = pf.guess[i] + cosmo.distmod(zcontrolCC).value
                        self.pardict[par]['prior_mean'] = pf.prior[i] + cosmo.distmod(zcontrolCC).value
                        self.pardict[par]['idx'] = idx + np.arange(len(zcontrolCC))
                        self.pardict[par]['bounds'] = (pf.lbound[i] + cosmo.distmod(zcontrol).value,
                                                       pf.ubound[i] + cosmo.distmod(zcontrol).value)

                        if pf.use[i]: idx += len(zcontrolCC)
                elif pf.bins[i]:
                    if pf.bins[i] == 1:
                        self.pardict[par]['guess'] = np.zeros(len(zcontrol)) + pf.guess[i]
                        self.pardict[par]['prior_mean'] = np.zeros(len(zcontrol)) + pf.prior[i]
                        self.pardict[par]['idx'] = idx + np.arange(len(zcontrol))
                        if pf.use[i]: idx += len(zcontrol)
                    else:
                        zcontrolCC = np.logspace(np.log10(min(zcontrol)),np.log10(max(zcontrol)),pf.bins[i])
                        self.pardict[par]['guess'] = np.zeros(len(zcontrolCC)) + pf.guess[i]
                        self.pardict[par]['prior_mean'] = np.zeros(len(zcontrolCC)) + pf.prior[i]
                        self.pardict[par]['idx'] = idx + np.arange(len(zcontrolCC))
                        if pf.use[i]: idx += len(zcontrolCC)
                elif pf.zpoly[i]:
                    self.pardict[par]['guess'] = np.append(pf.guess[i],np.array([0.]*int(pf.zpoly[i])))
                    self.pardict[par]['prior_mean'] = np.append(pf.prior[i],np.array([0.]*int(pf.zpoly[i])))
                    self.pardict[par]['idx'] = idx + np.arange(int(pf.zpoly[i])+1)
                    if pf.use[i]: idx += int(pf.zpoly[i])+1
                elif pf.use[i]: idx += 1

                if pf.fixed[i]:
                    self.pardict[par]['prior_std'] = 1e-5
#.........這裏部分代碼省略.........

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import distmod [as 別名]
def zmodel(x,zcontrol,zHD,pardict,corr=True):
    if not corr: from astropy.cosmology import Planck13 as cosmo
        
    muAmodel = np.zeros(len(zHD))
    if pardict['popBmean']['bins'] or pardict['popBmean']['zpoly']:
        muBmodel = np.zeros(len(zHD))
    else: muBmodel = None
    if pardict['popBstd']['bins'] or pardict['popBstd']['zpoly']:
        sigBmodel = np.zeros(len(zHD))
    else: sigBmodel = None
    if pardict['popB2mean']['bins'] or pardict['popB2mean']['zpoly']:
        muB2model = np.zeros(len(zHD))
    else: muB2model = None
    if pardict['popB2std']['bins'] or pardict['popB2std']['zpoly']:
        sigB2model = np.zeros(len(zHD))
    else: sigB2model = None
    if pardict['skewB']['bins'] or pardict['skewB']['zpoly']:
        skewBmodel = np.zeros(len(zHD))
    else: skewBmodel = None

    # Ia redshift/distance model
    for zb,zb1,i in zip(zcontrol[:-1],zcontrol[1:],range(len(zcontrol))):
        mua,mua1 = x[pardict['popAmean']['idx'][i]],x[pardict['popAmean']['idx'][i+1]]
        cols = np.where((zHD >= zb) & (zHD < zb1))[0]
        if not corr: cosmod = cosmo.distmod(zHD[cols]).value; cosbin = cosmo.distmod(zb).value
        alpha = np.log10(zHD[cols]/zb)/np.log10(zb1/zb)
        if corr:
            muAmodel[cols] = (1-alpha)*mua + alpha*mua1
        else:
            muAmodel[cols] = mua + cosmod - cosbin

    # CC redshift/distance model - need same # of bins for everything CC-related ATM
    if pardict['popBmean']['use'] and pardict['popBmean']['bins']:
        zcontrolCC = np.logspace(np.log10(min(zcontrol)),np.log10(max(zcontrol)),len(pardict['popBmean']['idx']))
        for zb,zb1,i in zip(zcontrolCC[:-1],zcontrolCC[1:],range(len(zcontrol))):
            cols = np.where((zHD >= zb) & (zHD < zb1))[0]
            if not corr: cosmod = cosmo.distmod(zHD[cols]).value; cosbin = cosmo.distmod(zb).value
            alpha = np.log10(zHD[cols]/zb)/np.log10(zb1/zb)
            if pardict['popBmean']['use'] and pardict['popBmean']['bins']:
                mub_1,mub1_1 = x[pardict['popBmean']['idx'][i]],x[pardict['popBmean']['idx'][i+1]]
                if corr:
                    muBmodel[cols] = (1-alpha)*mub_1 + alpha*mub1_1
                else:
                    muBmodel[cols] = mub_1 + cosmod - cosbin
            if pardict['popBstd']['use'] and pardict['popBstd']['bins']:
                sigb_1,sigb1_1 = x[pardict['popBstd']['idx'][i]],x[pardict['popBstd']['idx'][i+1]]
                if corr:
                    sigBmodel[cols] = (1-alpha)*sigb_1 + alpha*sigb1_1
                else:
                    sigBmodel[cols] = sigb_1
            if pardict['skewB']['use'] and pardict['skewB']['bins']:
                skewb,skewb1 = x[pardict['skewB']['idx'][i]],x[pardict['skewB']['idx'][i+1]]
                if corr:
                    skewBmodel[cols] = (1-alpha)*skewb + alpha*skewb1
                else:
                    skewBmodel[cols] = skewb

    # second gaussian - allowing for different # of bins
    if pardict['popB2mean']['use'] and pardict['popB2mean']['bins']:
        zcontrolCC = np.logspace(np.log10(min(zcontrol)),np.log10(max(zcontrol)),len(pardict['popB2mean']['idx']))
        for zb,zb1,i in zip(zcontrolCC[:-1],zcontrolCC[1:],range(len(zcontrol))):
            cols = np.where((zHD >= zb) & (zHD < zb1))[0]
            if not corr: cosmod = cosmo.distmod(zHD[cols]).value; cosbin = cosmo.distmod(zb).value
            alpha = np.log10(zHD[cols]/zb)/np.log10(zb1/zb)
            if pardict['popB2mean']['use'] and pardict['popB2mean']['bins']:
                mub_2,mub1_2 = x[pardict['popB2mean']['idx'][i]],x[pardict['popB2mean']['idx'][i+1]]
                if corr:
                    muB2model[cols] = (1-alpha)*mub_2 + alpha*mub1_2
                else:
                    muB2model[cols] = mub_2 + cosmod - cosbin
            if pardict['popB2std']['use'] and pardict['popB2std']['bins']:
                sigb2,sigb1_2 = x[pardict['popB2std']['idx'][i]],x[pardict['popB2std']['idx'][i+1]]
                if corr:
                    sigB2model[cols] = (1-alpha)*sigb2 + alpha*sigb1_2
                else:
                    sigB2model[cols] = sigb2


    if pardict['skewB']['use'] and not pardict['skewB']['bins'] and not pardict['skewB']['zpoly']:
        skewBmodel = np.zeros(len(zHD)) + x[pardict['skewB']['idx']]
    if pardict['popBmean']['use'] and not pardict['popBmean']['bins'] and not pardict['popBmean']['zpoly']:
        muBmodel = muAmodel + x[pardict['popBmean']['idx']]
    if pardict['popBstd']['use'] and not pardict['popBstd']['bins'] and not pardict['popBstd']['zpoly']:
        sigBmodel = x[pardict['popBstd']['idx']]
    if pardict['popBmean']['use'] and pardict['popBmean']['zpoly']:
        muBmodel = 1*muAmodel
        sigBmodel = 0
        for i,j in zip(pardict['popBmean']['idx'],
                       range(len(pardict['popBmean']['idx'])-1)):
            muBmodel += x[pardict['popBmean']['idx'][j]]*zHD**j/(1+pardict['popBmean']['idx'][-1]*zHD)
    if pardict['popBstd']['use'] and pardict['popBstd']['zpoly']:
        for i,j in zip(pardict['popBstd']['idx'],
                       range(len(pardict['popBstd']['idx']))):
            sigBmodel += x[pardict['popBstd']['idx'][j]]*zHD**j/(1+pardict['popBstd']['idx'][-1]*zHD)

    if pardict['popB2mean']['use'] and not pardict['popB2mean']['bins'] and not pardict['popB2mean']['zpoly']:
        muB2model = muAmodel + x[pardict['popB2mean']['idx']]
    if pardict['popB2std']['use'] and not pardict['popB2std']['bins'] and not pardict['popB2std']['zpoly']:
        sigB2model = x[pardict['popB2std']['idx']]
    if pardict['popB2mean']['use'] and pardict['popB2mean']['zpoly']:
#.........這裏部分代碼省略.........

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import distmod [as 別名]
    def mcmc(self,inp,zcontrol,guess):
        from scipy.optimize import minimize,basinhopping
        import emcee
        if not inp.__dict__.has_key('PL'):
            inp.PL = 0

        # minimize, not maximize
        if self.pardict['popB2mean']['use']:
            lnlikefunc = lambda *args: -threegausslike(*args)
        elif self.pardict['skewB']['use']:
            lnlikefunc = lambda *args: -twogausslike_skew(*args)
        else:
            lnlikefunc = lambda *args: -twogausslike(*args)

        inp.mu,inp.muerr = salt2mu(x1=inp.x1,x1err=inp.x1ERR,c=inp.c,cerr=inp.cERR,mb=inp.mB,mberr=inp.mBERR,
                                   cov_x1_c=inp.COV_x1_c,cov_x1_x0=inp.COV_x1_x0,cov_c_x0=inp.COV_c_x0,
                                   alpha=self.options.salt2alpha,
                                   beta=self.options.salt2beta,
                                   x0=inp.x0,z=inp.zHD)
        inp.residerr = inp.muerr[:]
        inp.resid = inp.mu - cosmo.distmod(inp.zHD).value

                        
        bounds,usebounds = (),False
        for i in xrange(len(guess)): 
            key = getpar(i,self.pardict)
            if (not self.pardict[key]['addcosmo'] and self.pardict[key]['bounds'][0] != self.pardict[key]['bounds'][1]) or \
                    (self.pardict[key]['addcosmo'] and self.pardict[key]['bounds'][0][0] != self.pardict[key]['bounds'][1][0]): 
                if hasattr(self.pardict[key]['idx'],"__len__"):
                    if i in self.pardict[key]['idx']:
                        if self.pardict[key]['addcosmo']:
                            bounds += ((self.pardict[key]['bounds'][0][self.pardict[key]['idx'] == i][0],
                                        self.pardict[key]['bounds'][1][self.pardict[key]['idx'] == i][0]),)
                        else:
                            bounds += (self.pardict[key]['bounds'],)
                        usebounds = True
                    else:
                        bounds += ((None,None),)
                else:
                    if i == self.pardict[key]['idx']:
                        bounds += (self.pardict[key]['bounds'],)
                        usebounds = True
                    else:
                        bounds += ((None,None),)
            else:
                bounds += ((None,None),)
        if not usebounds: bounds = None

        if not self.options.ntemps:
            if self.options.miniter <= 1:
                md = minimize(lnlikefunc,guess,
                              args=(inp,zcontrol,self.pardict['scaleA']['use'],self.pardict,self.options.debug),
                              bounds=bounds,method=self.options.minmethod,options={'maxiter':10000,'maxfev':10000})
            else:
                md = basinhopping(lnlikefunc,guess,
                                  minimizer_kwargs = {'args':(inp,zcontrol,self.pardict['scaleA']['use'],
                                                              self.pardict,self.options.debug),
                                                      'method':self.options.minmethod,
                                                      'bounds':bounds},niter=self.options.miniter)

            if md.message != 'Optimization terminated successfully.' and \
                    md.message != 'requested number of basinhopping iterations completed successfully':
                print(md.message)
                if self.options.forceminsuccess:
                    raise RuntimeError('Error : Minimization Failed!!!')
                else:
                    print("""Warning : Minimization Failed!!!
Try some different initial guesses, or let the MCMC try and take care of it""")

        # add in the random steps
        if not self.options.ntemps:
            ndim, nwalkers = len(md["x"]), int(self.options.nwalkers)
            mcstep = np.array([])
            for i in range(len(md["x"])):
                mcstep = np.append(mcstep,getparval(i,self.pardict,'mcstep'))
            pos = [md["x"] + mcstep*np.random.randn(ndim) for i in range(nwalkers)]
        else:
            ndim, nwalkers = len(guess), int(self.options.nwalkers)
            mcstep = np.array([])
            for i in range(len(guess)):
                mcstep = np.append(mcstep,getparval(i,self.pardict,'mcstep'))
            pos = [np.array(guess) + mcstep*np.random.randn(ndim) for i in range(nwalkers)]


        if not self.options.ntemps:
            sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob,
                                            args=(inp,zcontrol,self.pardict,self.options.debug),
                                            threads=int(self.options.nthreads))
            pos, prob, state = sampler.run_mcmc(pos, self.options.ninit)
            sampler.reset()
            sampler.run_mcmc(pos, self.options.nsteps, thin=1)
        else:
            if self.pardict['popB2mean']['use']:
                lnlikefunc = threegausslike
            elif self.pardict['skewB']['use']:
                lnlikefunc = twogausslike_skew
            else:
                lnlikefunc = twogausslike

            sampler = emcee.PTSampler(self.options.ntemps,nwalkers, ndim, lnlikefunc,
#.........這裏部分代碼省略.........