Python maxentropy.logsumexp函数代码示例

def calc_avg_rcmks(parser):
    options,args= parser.parse_args()
    njks= 101
    nmks= 101
    jks= numpy.linspace(0.5,0.8,njks)
    mks= numpy.linspace(-0.5,-3.,nmks)
    if options.basti:
        zs= numpy.array([0.004,0.008,0.01,0.0198,0.03,0.04])
        zsolar= 0.019
    elif options.parsec:
        zs= numpy.arange(0.0005,0.06005,0.0005)
#        zs= numpy.array([0.01,0.02])
        zsolar= 0.019
    else:
        zs= numpy.arange(0.0005,0.03005,0.0005)
#        zs= numpy.array([0.01,0.02])
        zsolar= 0.019
    if not os.path.exists(options.outfilename):
        logpz= localzdist(zs,zsolar=zsolar)
        logmkp= numpy.zeros((len(zs),njks,nmks))
        logp= numpy.zeros((len(zs),njks,nmks))      
        funcargs= (zs,options,njks,jks,nmks,mks,logpz)
        multOut= multi.parallel_map((lambda x: indiv_calc(x,
                                                          *funcargs)),
                                    range(len(zs)),
                                    numcores=numpy.amin([64,len(zs),
                                                         multiprocessing.cpu_count()]))
        for ii in range(len(zs)):
            logmkp[ii,:,:]= multOut[ii][0,:,:]
            logp[ii,:,:]= multOut[ii][1,:,:]
        save_pickles(options.outfilename,logmkp,logp)
    else:
        savefile= open(options.outfilename,'rb')
        logmkp= pickle.load(savefile)
        logp= pickle.load(savefile)
        savefile.close()
    indx= numpy.isnan(logp)
    logp[indx]= -numpy.finfo(numpy.dtype(numpy.float64)).max
    logmkp[indx]= -numpy.finfo(numpy.dtype(numpy.float64)).max
    #Average the peak, so calculate the peak
    for ii in range(len(zs)):
        for jj in range(njks):
            maxmkindx= numpy.argmax(logp[ii,jj,:])
            totlogp= maxentropy.logsumexp(logp[ii,jj,:])
            logmkp[ii,jj,:]= logmkp[ii,jj,maxmkindx]-logp[ii,jj,maxmkindx]+totlogp
            logp[ii,jj,:]= totlogp
    avgmk= numpy.exp(maxentropy.logsumexp(logmkp.flatten())\
                         -maxentropy.logsumexp(logp.flatten()))
    solindx= numpy.argmin(numpy.fabs(zs-0.017))
    avgmksolar= numpy.exp(maxentropy.logsumexp(logmkp[solindx,:,:].flatten())\
                              -maxentropy.logsumexp(logp[solindx,:,:].flatten()))
    print "Average mk: %f" % (-avgmk)
    print "Average mk if solar: %f" % (-avgmksolar)
    return -avgmk

def calc_model(params,options,data,logpiso,logpisodwarf,df,nlocs,locations,iso):
    avg_plate_model= numpy.zeros(nlocs)
    for ii in range(nlocs):
        #Calculate vlos | los
        indx= (data['LOCATION'] == locations[ii])
        thesedata= data[indx]
        thislogpiso= logpiso[indx,:]
        if options.dwarf:
            thislogpisodwarf= logpisodwarf[indx,:]
        else:
            thislogpisodwarf= None
        vlos= numpy.linspace(-200.,200.,options.nvlos)
        pvlos= numpy.zeros(options.nvlos)
        if not options.multi is None:
            pvlos= multi.parallel_map((lambda x: pvlosplate(params,vlos[x],
                                                            thesedata,
                                                            df,options,
                                                            thislogpiso,
                                                            thislogpisodwarf,iso)),
                                      range(options.nvlos),
                                      numcores=numpy.amin([len(vlos),multiprocessing.cpu_count(),options.multi]))
        else:
            for jj in range(options.nvlos):
                print jj
                pvlos[jj]= pvlosplate(params,vlos[jj],thesedata,df,options,
                                      thislogpiso,thislogpisodwarf)
        pvlos-= logsumexp(pvlos)
        pvlos= numpy.exp(pvlos)
        #Calculate mean and velocity dispersion
        avg_plate_model[ii]= numpy.sum(vlos*pvlos)
    return avg_plate_model

 def call_polymorphism(self, obs, post):
     """Get the polymorphism probability.
     This is the posterior probability that the strain is homozygous
     for the non-reference base with the highest count at this position.
     @param obs: one ref base count and three non-ref base counts
     @param post: the posterior hidden state probabilities
     @return: the polymorphism probability
     """
     # unpack the posterior state distribution
     p_recent, p_ancient, p_garbage, p_misaligned = post
     # get the prior probability of polymorphism conditional on state
     p_recent_AA = self.states[0].get_posterior_distribution(obs)[2]
     p_ancient_AA = self.states[1].get_posterior_distribution(obs)[2]
     # compute the posterior probability of a polymorphism
     posterior_polymorphism = 0
     posterior_polymorphism += p_recent * p_recent_AA
     posterior_polymorphism += p_ancient * p_ancient_AA
     # Given that a polymorphism occurred,
     # get the probability distribution over the
     # three non-reference nucleotides.
     r = self.seqerr
     log_Pr = math.log(r/4.0)
     log_PA = math.log(1 - 3*r/4.0)
     logs = [
             obs[1]*log_PA + obs[2]*log_Pr + obs[3]*log_Pr,
             obs[1]*log_Pr + obs[2]*log_PA + obs[3]*log_Pr,
             obs[1]*log_Pr + obs[2]*log_Pr + obs[3]*log_PA]
     condmaxpost = math.exp(max(logs) - logsumexp(logs))
     # get the posterior probability distribution
     maxpost = posterior_polymorphism * condmaxpost
     return maxpost

def run(*args):
	dprintn(8, "# Generating data")
	
	global hypotheses
	RANK = str(MPI.COMM_WORLD.Get_rank())
	
	data_size = args[0]
	
	p_representation = defaultdict(int) # how often do you get the right representation
	p_response = defaultdict(int) # how often do you get the right response?
	p_representation_literal = defaultdict(int) # how often do you get the right representation
	p_response_literal = defaultdict(int)  # how often do you get the right response?
	p_representation_presup = defaultdict(int) # how often do you get the right representation
	p_response_presup = defaultdict(int) # how often do you get the right response?
	
	dprintn(8, "# Generating data")
	data = generate_data(data_size)
	
	# recompute these
	dprintn(8, "# Computing posterior")
	#[ x.unclear_functions() for x in hypotheses ]
	[ x.compute_posterior(data) for x in hypotheses ]
	
	# normalize the posterior in fs
	dprintn(8, "# Computing normalizer")
	Z = logsumexp([x.posterior_score for x in hypotheses])
	
	# and output the top hypotheses
	qq = FiniteBestSet(max=True, N=25)
	for h in hypotheses: qq.push(h, h.posterior_score) # get the tops
	for i, h in enumerate(qq.get_sorted()):
		for w in h.all_words():
			fprintn(8, data_size, i, w, h.posterior_score, q(h.lex[w]), f=options.OUT_PATH+"-hypotheses."+RANK+".txt")
	
	# and compute the probability of being correct
	dprintn(8, "# Computing correct probability")
	for h in hypotheses:
		hstr = str(h)
		#print data_size, len(data), exp(h.posterior_score), correct[ str(h)+":"+w ]
		for w in words:
			p = exp(h.posterior_score - Z)
			key = w + ":" + hstr
			
			p_representation[w] += p * (agree_pct[key] == 1.)
			p_representation_presup[w]  += p * (agree_pct_presup[key] == 1.) # if we always agree with the target, then we count as the right rep.
			p_representation_literal[w] += p * (agree_pct_literal[key] == 1.)
			
			# and just how often does the hypothesis agree?
			p_response[w] += p * agree_pct[key]
			p_response_presup[w]  += p * agree_pct_presup[key]
			p_response_literal[w] += p * agree_pct_literal[key]
			
	dprintn(8, "# Outputting")
	

	for w in words:
		fprintn(10, rank, q(w), data_size, p_representation[w], p_representation_presup[w], p_representation_literal[w], p_response[w], p_response_presup[w], p_response_literal[w], f=options.OUT_PATH+"-stats."+RANK+".txt")
	
	return 0

def testErrs(options,args):
    ndfehs, ndafes= 201,201
    dfehs= numpy.linspace(0.01,0.4,ndfehs)
    dafes= numpy.linspace(0.01,0.3,ndafes)
    if os.path.exists(args[0]):
        savefile= open(args[0],'rb')
        loglike= pickle.load(savefile)
        ii= pickle.load(savefile)
        jj= pickle.load(savefile)
        savefile.close()
    else:
        loglike= numpy.zeros((ndfehs,ndafes))
        ii, jj= 0, 0
    while ii < ndfehs:
        while jj < ndafes:
            sys.stdout.write('\r'+"Working on %i / %i" %(ii*ndafes+jj+1,ndafes*ndfehs))
            sys.stdout.flush()
            loglike[ii,jj]= errsLogLike(dfehs[ii],dafes[jj],options)
            jj+= 1
        ii+= 1
        jj= 0
        save_pickles(args[0],loglike,ii,jj)
    save_pickles(args[0],loglike,ii,jj)
    sys.stdout.write('\r'+_ERASESTR+'\r')
    sys.stdout.flush()
    if options.prior:
        prior= numpy.zeros((ndfehs,ndafes))
        for ii in range(ndfehs):
            prior[ii,:]= -0.5*(dafes-0.1)**2./0.1**2.-0.5*(dfehs[ii]-0.2)**2./0.1**2.
        loglike+= prior
    loglike-= maxentropy.logsumexp(loglike)
    loglike= numpy.exp(loglike)
    loglike/= numpy.sum(loglike)*(dfehs[1]-dfehs[0])*(dafes[1]-dafes[0])
    #Plot
    bovy_plot.bovy_print()
    bovy_plot.bovy_dens2d(loglike.T,origin='lower',
                          cmap='gist_yarg',
                          xlabel=r'\delta_{[\mathrm{Fe/H}]}',
                          ylabel=r'\delta_{[\alpha/\mathrm{Fe}]}',
                          xrange=[dfehs[0],dfehs[-1]],
                          yrange=[dafes[0],dafes[-1]],
                          contours=True,
                          cntrmass=True,
                          onedhists=True,
                          levels= special.erf(0.5*numpy.arange(1,4)))
    if options.prior:
        bovy_plot.bovy_text(r'$\mathrm{with\ Gaussian\ prior:}$'+
                            '\n'+r'$\delta_{[\mathrm{Fe/H}]}= 0.2 \pm 0.1$'
                            +'\n'+r'$\delta_{[\alpha/\mathrm{Fe}]}= 0.1 \pm 0.1$',
top_right=True)
    bovy_plot.bovy_end_print(options.plotfile)

    def neg_log_likelihood(theta_sparse, hb = None):
        if not hb is None:
            h, b = hb
        else:
            h, b = dp(theta_sparse)
        
        log_kappa = logsumexp(h[0] + b[1])

        nll = log_kappa
        nll -= h[0][0]
        for k in range(1, params['M']):
            nll -= h[k][0,0]
        for ind in theta_sparse:
            nll += params['lambda'] * np.abs(theta_sparse[ind])
        return nll

def pvlosplate(params,vhelio,data,df,options,logpiso,logpisodwarf,iso):
    """
    NAME:
       pvlosplate
    PURPOSE:
       calculate the vlos probability for a given location
    INPUT:
       params - parameters of the model
       vhelio - heliocentric los velocity to evaluate
       data - data array for this location
       df - df object(s) (?)
       options - options
       logpiso, logpisodwarf - precalculated isochrones
    OUTPUT:
       log of the probability
    HISTORY:
       2012-02-20 - Written - Bovy (IAS)
    """
    #Output is sum over data l,b,jk,h
    l= data['GLON']*_DEGTORAD
    b= data['GLAT']*_DEGTORAD
    sinl= numpy.sin(l)
    cosl= numpy.cos(l)
    sinb= numpy.sin(b)
    cosb= numpy.cos(b)
    jk= data['J0MAG']-data['K0MAG']
    try:
        jk[(jk < 0.5)]= 0.5 #BOVY: FIX THIS HACK BY EMAILING GAIL
    except TypeError:
        pass #HACK
    h= data['H0MAG']
    options.multi= 1 #To avoid conflict
    out= -mloglike(params,numpy.zeros(len(data))+vhelio,
                   l,
                   b,
                   jk,
                   h,
                   df,options,
                   sinl,
                   cosl,
                   cosb,
                   sinb,
                   logpiso,
                   logpisodwarf,True,None,iso,data['FEH']) #None iso for now
    #indx= (out >= -0.1)*(out <= 0.1)
    #print out[indx], jk[indx], h[indx]
    return logsumexp(out)

    def bound(self, corpus, gamma=None, subsample_ratio=1.0):
        """
        Estimate the variational bound of documents from `corpus`:
        E_q[log p(corpus)] - E_q[log q(corpus)]

        `gamma` are the variational parameters on topic weights for each `corpus`
        document (=2d matrix=what comes out of `inference()`).
        If not supplied, will be inferred from the model.

        """
        score = 0.0
        _lambda = self.state.get_lambda()
        Elogbeta = dirichlet_expectation(_lambda)

        for d, doc in enumerate(corpus):  # stream the input doc-by-doc, in case it's too large to fit in RAM
            if d % self.chunksize == 0:
                logger.debug("bound: at document #%i", d)
            if gamma is None:
                gammad, _ = self.inference([doc])
            else:
                gammad = gamma[d]
            Elogthetad = dirichlet_expectation(gammad)

            # E[log p(doc | theta, beta)]
            score += np.sum(cnt * logsumexp(Elogthetad + Elogbeta[:, int(id)]) for id, cnt in doc)

            # E[log p(theta | alpha) - log q(theta | gamma)]; assumes alpha is a vector
            score += np.sum((self.alpha - gammad) * Elogthetad)
            score += np.sum(gammaln(gammad) - gammaln(self.alpha))
            score += gammaln(np.sum(self.alpha)) - gammaln(np.sum(gammad))

        # Compensate likelihood for when `corpus` above is only a sample of the whole corpus. This ensures
        # that the likelihood is always rougly on the same scale.
        score *= subsample_ratio

        # E[log p(beta | eta) - log q (beta | lambda)]; assumes eta is a scalar
        score += np.sum((self.eta - _lambda) * Elogbeta)
        score += np.sum(gammaln(_lambda) - gammaln(self.eta))

        if np.ndim(self.eta) == 0:
            sum_eta = self.eta * self.num_terms
        else:
            sum_eta = np.sum(self.eta)

        score += np.sum(gammaln(sum_eta) - gammaln(np.sum(_lambda, 1)))

        return score

 def _parse_hz_dict_indiv(self,hz):
     htype= hz.get('type','exp')
     if htype == 'exp':
         zd= hz.get('h',0.0375)
         th= lambda z, tzd=zd: 1./2./tzd*numpy.exp(-numpy.fabs(z)/tzd)
         tH= lambda z, tzd= zd: (numpy.exp(-numpy.fabs(z)/tzd)-1.
                                 +numpy.fabs(z)/tzd)*tzd/2.
         tdH= lambda z, tzd= zd: 0.5*numpy.sign(z)\
             *(1.-numpy.exp(-numpy.fabs(z)/tzd))
     elif htype == 'sech2':
         zd= hz.get('h',0.0375)
         th= lambda z, tzd=zd: 1./numpy.cosh(z/2./tzd)**2./4./tzd
         # Avoid overflow in cosh
         tH= lambda z, tzd= zd: \
             tzd*(logsumexp(numpy.array([z/2./tzd,-z/2./tzd]),axis=0)\
                      -numpy.log(2.))
         tdH= lambda z, tzd= zd: numpy.tanh(z/2./tzd)/2.
     return (th,tH,tdH)

def _eval_sumgaussians(x,xamp,xmean,xcovar):
    """x array [ndata,ndim], return log"""
    ndata= x.shape[0]
    da= x.shape[1]
    out= numpy.zeros(ndata)
    ngauss= len(xamp)
    loglike= numpy.zeros(ngauss)
    for ii in range(ndata):
        for kk in range(ngauss):
            if xamp[kk] == 0.:
                loglike[kk]= numpy.finfo(numpy.dtype(numpy.float64)).min
                continue
            tinv= linalg.inv(xcovar[kk,:,:])
            delta= x[ii,:]-xmean[kk,:]
            loglike[kk]= numpy.log(xamp[kk])+0.5*numpy.log(linalg.det(tinv))\
                -0.5*numpy.dot(delta,numpy.dot(tinv,delta))+\
                da*_SQRTTWOPI
        out[ii]= maxentropy.logsumexp(loglike)
    return out

 def bound(self, corpus, gamma=None, subsample_ratio=1.0):
     score = 0.0
     _lambda = self.state.get_lambda()
     Elogbeta = dirichlet_expectation(_lambda)
     for d, doc in enumerate(corpus):
         if gamma is None:
             gammad, _ = self.inference([doc])
         else:
             gammad = gamma[d]
         Elogthetad = dirichlet_expectation(gammad)
         score += numpy.sum(cnt * logsumexp(Elogthetad + Elogbeta[:, id]) for id, cnt in doc)
         score += numpy.sum((self.alpha - gammad) * Elogthetad)
         score += numpy.sum(gammaln(gammad) - gammaln(self.alpha))
         score += gammaln(numpy.sum(self.alpha)) - gammaln(numpy.sum(gammad))
     score *= subsample_ratio
     score += numpy.sum((self.eta - _lambda) * Elogbeta)
     score += numpy.sum(gammaln(_lambda) - gammaln(self.eta))
     score += numpy.sum(gammaln(self.eta * self.num_terms) - gammaln(numpy.sum(_lambda, 1)))
     return score

    def inference(self, doc):
        """
        Perform inference on a single document.
        
        Return 3-tuple of `(likelihood of this document, word-topic distribution
        phi, expected word counts gamma (~topic distribution))`.
        
        A document is simply a bag-of-words collection which supports len() and 
        iteration over (wordIndex, wordCount) 2-tuples.
        
        The model itself is not affected in any way (this function is read-only aka 
        const).
        """
        # init help structures
        totalWords = sum(wordCount for _, wordCount in doc)
        gamma = numpy.zeros(self.numTopics) + self.alpha + 1.0 * totalWords / self.numTopics
        phi = numpy.zeros(shape = (len(doc), self.numTopics)) + 1.0 / self.numTopics
        likelihood = likelihoodOld = converged = numpy.NAN
        
        # variational estimate
        for i in xrange(self.VAR_MAX_ITER):
#            logger.debug("inference step #%s, converged=%s, likelihood=%s, likelikelihoodOld=%s" % 
#                          (i, converged, likelihood, likelihoodOld))
            
            if numpy.isfinite(converged) and converged <= self.VAR_CONVERGED:
                logger.debug("document converged in %i iterations" % i)
                break
            
            for n, (wordIndex, wordCount) in enumerate(doc):
                # compute phi vars, in log space, to prevent numerical nastiness
                tmp = digamma(gamma) + self.logProbW[:, wordIndex] # vector operation

                # convert phi and update gamma
                newPhi = numpy.exp(tmp - logsumexp(tmp))
                gamma += wordCount * (newPhi - phi[n])
                phi[n] = newPhi
            
            likelihood = self.computeLikelihood(doc, phi, gamma)
            assert numpy.isfinite(likelihood)
            converged = numpy.divide(likelihoodOld - likelihood, likelihoodOld)
            likelihoodOld = likelihood
        return likelihood, phi, gamma

    def bound(self, corpus, gamma=None):
        """
        Estimate the variational bound of documents from `corpus`.

        `gamma` are the variational parameters on topic weights (one for each
        document in `corpus`). If not supplied, will be automatically inferred
        from the model.
        """
        score = 0.0
        Elogbeta = numpy.log(self.expElogbeta)

        for d, doc in enumerate(corpus):
            if d % self.chunks == 0:
                logger.info("PROGRESS: at document #%i" % d)
            if gamma is None:
                gammad, _ = self.inference([doc])
            else:
                gammad = gamma[d, :]
            Elogthetad = dirichlet_expectation(gammad)
            expElogthetad = numpy.exp(Elogthetad)
            ids = [id for id, _ in doc]
            cts = numpy.array([cnt for _, cnt in doc])
            phinorm = numpy.zeros(len(ids))
            for i in xrange(len(ids)):
                phinorm[i] = logsumexp(Elogthetad + Elogbeta[:, ids[i]])

            # E[log p(docs | theta, beta)]
            score += numpy.sum(cts * phinorm)

            # E[log p(theta | alpha) - log q(theta | gamma)]
            score += numpy.sum((self.alpha - gammad) * Elogthetad)
            score += numpy.sum(gammaln(gammad) - gammaln(self.alpha))
            score += gammaln(self.alpha * self.numTopics) - gammaln(numpy.sum(gammad))

        # E[log p(beta | eta) - log q (beta | lambda)]
        score += numpy.sum((self.eta - self._lambda) * Elogbeta)
        score += numpy.sum(gammaln(self._lambda) - gammaln(self.eta))
        score += numpy.sum(gammaln(self.eta * self.numTerms) -
                              gammaln(numpy.sum(self._lambda, 1)))

        return score

    def dp(theta_sparse):
        theta = theta_dense(theta_sparse)

        h = [None] * params['M']
        h[0] = np.empty(n_w[0])
        for w in range(n_w[0]):
            h[0][w] = np.sum(theta * hits_pre[0][w])
        for k in range(1, params['M']):
            h[k] = np.empty((n_w[k-1], n_w[k]))
            for w_prev in range(n_w[k-1]):
                for w in range(n_w[k]):
                    h[k][w_prev,w] = np.sum(theta * hits_pre[k][w_prev,w])

        b = [None] * (params['M']+1)
        b[params['M']] = np.zeros(n_w[params['M']-1])
        for k in range(params['M']-1, 0, -1):
            b[k] = np.empty(n_w[k-1])
            for w_prev in range(n_w[k-1]):
                b[k][w_prev] = logsumexp(h[k][w_prev] + b[k+1])

        return h, b

def _eval_gauss_grid(x,y,xamp,xmean,xcovar):
    nx= len(x)
    ny= len(y)
    out= numpy.zeros((nx,ny))
    ngauss= len(xamp)
    dim= xmean.shape[1]
    loglike= numpy.zeros(ngauss)
    for ii in range(nx):
        for jj in range(ny):
            a= numpy.array([x[ii],y[jj]])
            for kk in range(ngauss):
                if xamp[kk] == 0.:
                    loglike[kk]= numpy.finfo(numpy.dtype(numpy.float64)).min
                    continue
                tinv= numpy.linalg.inv(xcovar[kk,:,:])
                delta= a-xmean[kk,:]
                loglike[kk]= numpy.log(xamp[kk])+0.5*numpy.log(numpy.linalg.det(tinv))\
                    -0.5*numpy.dot(delta,numpy.dot(tinv,delta))+\
                    dim*_SQRTTWOPI
            out[ii,jj]= logsumexp(loglike)
    return out