def compute_likelihood(document, model, phi, var_gamma):
    likelihood = 0
    digsum = 0
    var_gamma_sum = 0
    dig = [0 for x in range(model.num_topics)]

    for k in range(0, model.num_topics):
        dig[k] = digamma(var_gamma[k])
        var_gamma_sum = var_gamma[k] + var_gamma_sum

    digsum = digamma(var_gamma_sum)

    likelihood = math.lgamma(model.alpha * model.num_topics) \
                 - model.num_topics * math.lgamma(model.alpha) \
                 - (math.lgamma(var_gamma_sum))

    for k in range(0, model.num_topics):
        likelihood += ((model.alpha - 1) * (dig[k] - digsum)
                       + math.lgamma(var_gamma[k]) - (var_gamma[k] - 1)
                       * (dig[k] - digsum))

        for n in range(0, document.unique_word_count):
            if phi[n][k] > 0:
                likelihood += document.word_counts[n] * \
                              (phi[n][k] * ((dig[k] - digsum)
                                            - math.log(phi[n][k])
                                            + model.log_prob_w[k][document.words[n]]))

    return likelihood

 def log_likelihood(self, full=False):
     ll = (math.lgamma(self.alpha) - math.lgamma(self.alpha + self.total_customers)
             + sum(math.lgamma(c) for tables in self.tables.itervalues() for c in tables)
             + self.ntables * math.log(self.alpha))
     if full:
         ll += self.base.log_likelihood(full=True) + self.prior.log_likelihood()
     return ll

def LogCombinations(x,y):
    u"""Calculates the logarithm of a binomial coefficient.
    This avoids overflows. Implemented with gamma functions for efficiency"""
    result=lgamma(x+1)
    result-=lgamma(y+1)
    result-=lgamma(x-y+1)
    return result

 def log_likelihood(self, full=False):
     ll = (math.lgamma(self.K * self.alpha) - math.lgamma(self.K * self.alpha + self.N)
             + sum(math.lgamma(self.alpha + self.count[k]) for k in xrange(self.K))
             - self.K * math.lgamma(self.alpha))
     if full:
         ll += self.prior.log_likelihood()
     return ll

    def UpdateKappaSigmaSq(self,it):
        for ii in xrange(self.T-1):
            new_kappa_sigma_sq = self.kappa_sigma_sq[ii]+(2*np.ceil(2*np.random.random())-3)*(np.random.geometric(1.0/(1+np.exp(self.log_kappa_sigma_sqq[ii])))-1)

            if new_kappa_sigma_sq <0:
                accept = 0
            else:
                lam1 = 1.0*self.lambda_sigma + self.kappa_sigma_sq[ii]
                gam1 = 1.0*self.lambda_sigma/self.mu_sigma + 1.0*self.rho_sigma/(1-self.rho_sigma)*self.lambda_sigma/self.mu_sigma
                loglike = lam1*np.log(gam1)-math.lgamma(lam1)+(lam1-1)*np.log(self.sigma_sq[ii+1])
                pnmean = self.sigma_sq[ii]*self.rho_sigma/(1-self.rho_sigma)*self.lambda_sigma/self.mu_sigma
                loglike = loglike + self.kappa_sigma_sq[ii]*np.log(pnmean)- math.lgamma(1.0*self.kappa_sigma_sq[ii]+1)

                lam1 = 1.0*self.lambda_sigma + new_kappa_sigma_sq
                gam1 = 1.0*self.lambda_sigma/self.mu_sigma + self.rho_sigma/(1-self.rho_sigma)*self.lambda_sigma/self.mu_sigma
                new_loglike = lam1*np.log(gam1)-math.lgamma(lam1)+(lam1-1)*np.log(self.sigma_sq[ii+1])
                pnmean = self.sigma_sq[ii]*self.rho_sigma/(1-self.rho_sigma)*self.lambda_sigma/self.mu_sigma
                new_loglike = new_loglike + new_kappa_sigma_sq*np.log(pnmean)-math.lgamma(1.0*new_kappa_sigma_sq+1)
                log_accept = new_loglike - loglike
                accept =1
                if np.isnan(log_accept) or np.isinf(log_accept):
                    accept = 0
                elif log_accept <0:
                    accept = np.exp(log_accept)

            self.kappa_lambda_sigma_accept = self.kappa_lambda_sigma_accept + accept
            self.kappa_sigma_sq_count = self.kappa_sigma_sq_count +1
            if np.random.random()<accept :
                self.kappa_sigma_sq[ii] = new_kappa_sigma_sq
            self.log_kappa_sigma_sqq[ii] = self.log_kappa_sigma_sqq[ii]+1.0/it**0.55*(accept-0.3)

 def sample_document(self, m):
     z = self.corpus[m]["state"]         # Step1: カウントを減らす
     if z > 0:
         self.topic_document_freq[z] -= 1
         self.topic_document_sum -= 1
         for v in self.corpus[m]["bag_of_words"]:
             self.topic_word_freq[z][v] -= 1
             self.topic_word_sum[z] -= 1
     n_d_v = defaultdict(float)          # Step2: 事後分布の計算
     n_d = 0.0
     for v in self.corpus[m]["bag_of_words"]:
         n_d_v[v] += 1.0
         n_d += 1.0
     p_z = defaultdict(lambda: 0.0)
     for z in xrange(1, self.K + 1):
         p_z[z] = math.log((self.topic_document_freq[z] + self.alpha) / (self.topic_document_sum + self.alpha*self.K))
         p_z[z] += (math.lgamma(self.topic_word_sum[z] + self.beta*self.V) - math.lgamma(self.topic_word_sum[z] + n_d + self.beta*self.V))
         for v in n_d_v.iterkeys():
             p_z[z] += (math.lgamma(self.topic_word_freq[z][v] + n_d_v[v] + self.beta) - math.lgamma(self.topic_word_freq[z][v] + self.beta))
     max_log = max(p_z.values())     # オーバーフロー対策
     for z in p_z:
         p_z[z] = math.exp(p_z[z] - max_log)
     new_z = self.sample_one(p_z)        # Step3: サンプル
     self.corpus[m]["state"] = new_z     # Step4: カウントを増やす
     self.topic_document_freq[new_z] += 1
     self.topic_document_sum += 1
     for v in self.corpus[m]["bag_of_words"]:
         self.topic_word_freq[new_z][v] += 1
         self.topic_word_sum[new_z] += 1

def incomplete_gamma(x, s):
    r"""
    This function computes the incomplete lower gamma function
    using the series expansion:

    .. math::

       \gamma(x, s) = x^s \Gamma(s)e^{-x}\sum^\infty_{k=0}
                    \frac{x^k}{\Gamma(s + k + 1)}

    This series will converge strongly because the Gamma
    function grows factorially.

    Because the Gamma function does grow so quickly, we can
    run into numerical stability issues. To solve this we carry
    out as much math as possible in the log domain to reduce
    numerical error. This function matches the results from
    scipy to numerical precision.
    """
    if x < 0:
        return 1
    if x > 1e3:
        return math.gamma(s)
    log_gamma_s = math.lgamma(s)
    log_x = log(x)
    value = 0
    for k in range(100):
        log_num = (k + s)*log_x + (-x) + log_gamma_s
        log_denom = math.lgamma(k + s + 1)
        value += math.exp(log_num - log_denom)
    return value

def calc_full(n, alphas):
    """ Calculate the log likelihood under DirMult distribution with alphas=avec, given data counts of nvec."""
    lg_sum_alphas = math.lgamma(alphas.sum())
    sum_lg_alphas = np.sum(scipy.special.gammaln(alphas))
    lg_sum_alphas_n = math.lgamma(alphas.sum() + n.sum())
    sum_lg_alphas_n = np.sum(scipy.special.gammaln(n+alphas))
    return lg_sum_alphas - sum_lg_alphas - lg_sum_alphas_n + sum_lg_alphas_n 

 def tdens(self, n, X):
     C = (1.0 + (X * X) / (n * 1.0))
     h = math.lgamma((n + 1.0) / 2.0) - math.lgamma(n / 2.0) 
     h = math.exp(h)
     h = h / math.sqrt(math.pi) / math.sqrt(n)
     Result = h * (C ** (-((n / 2.0) + (1.0 / 2.0))))
     return Result

def incompleteBetaFunction(x,a,b):
    lbeta = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) \
            + a * math.log(x) + b * math.log(1.0 - x)
    if (x < (a + 1)/(a + b + 2)):
        return math.exp(lbeta) * contFractionBeta(a,b,x)/a
    else:
        return 1 - math.exp(lbeta) * contFractionBeta(b,a,1.-x)/b

 def log_likelihood(self, full=False):
     ll = (math.lgamma(self.K * self.alpha) - math.lgamma(self.K * self.alpha + self.N)
             + sum(math.lgamma(self.alpha + self.count[k]) for k in self.count)
             - len(self.count) * math.lgamma(self.alpha)) # zero counts
     if full:
         ll += self.prior.log_likelihood()
     return ll

    def UpdateKappa(self, it):
        for ii in xrange(self.T-1):
            for jj in xrange(self.p):

                new_kappa = self.kappa[ii][jj]+(2*np.ceil(2*np.random.random())-3)*(np.random.geometric(1.0/(1+np.exp(self.log_kappa_q[ii][jj])))-1)

                if new_kappa < 0:
                    accept = 0
                else:
                    lam1 = self.lambda_[jj] + 1.0*self.kappa[ii][jj]
                    gam1 = self.lambda_[jj]/self.mu[jj] + self.delta[jj]
                    loglike = lam1*np.log(gam1) - math.lgamma(lam1)+(lam1-1)*np.log(self.psi[ii+1][jj])
                    pnmean = self.psi[ii][jj] * self.delta[jj]
                    loglike = loglike + 1.0*self.kappa[ii][jj]*np.log(pnmean) - math.lgamma(1.0*self.kappa[ii][jj]+1)

                    lam1 = self.lambda_[jj] + 1.0*new_kappa
                    gam1 = self.lambda_[jj]/self.mu[jj] + self.delta[jj]
                    new_loglike = lam1*np.log(gam1) - math.lgamma(lam1)+(lam1-1)*np.log(self.psi[ii+1][jj])
                    pnmean = self.psi[ii][jj]*self.delta[jj]
                    new_loglike = new_loglike + new_kappa*np.log(pnmean)-math.lgamma(1.0*new_kappa+1)
                    log_accept = new_loglike - loglike
                    accept =1
                    if np.isnan(log_accept) or np.isinf(log_accept):
                        accept =0
                    elif log_accept <0:
                        accept = np.exp(log_accept)

                self.kappa_accept = self.kappa_accept + accept
                self.kappa_count = self.kappa_count +1

                if np.random.random() < accept:
                    self.kappa[ii][jj] = new_kappa
                self.log_kappa_q[ii][jj] = self.log_kappa_q[ii][jj] + 1.0/it**0.55*(accept-0.3)

def theta_likelihood(theta, S, J):
    S += prior_s
    J += prior_j
    #If any of the values are 0 or negative return likelihood that will get rejected
    if theta <= 0 or S <= 0 or J <= 0:
        return 10000000
    else:
        return -(S * math.log(theta) + math.lgamma(theta) - math.lgamma(theta + J))

 def __compute_factor(self):
     self._factor = lgamma (self.community.J + 1)
     phi = table(self.community.abund)
     phi += [0] * int (max (self.community.abund) - len (phi))
     for spe in xrange (self.community.S):
         self._factor -= log (max (1, self.community.abund[spe]))
     for spe in xrange (int(max(self.community.abund))):
         self._factor -= lgamma (phi[spe] + 1)

 def logchoose(ni, ki):
     try:
         lgn1 = lgamma(ni + 1)
         lgk1 = lgamma(ki + 1)
         lgnk1 = lgamma(ni - ki + 1)
     except ValueError:
         raise ValueError
     return lgn1 - (lgnk1 + lgk1)

def multiTLogPDF(x,mu,Sigma,nu,p):

    part1 =   math.lgamma( 0.5 * (p + nu) )
    part2 = - math.lgamma( 0.5 * nu ) - 0.5 * p * np.log( nu ) - 0.5 * p * np.log( np.pi )
    part3 = - 0.5 * np.log( np.linalg.det(Sigma) )
    part4 = - 0.5 * ( nu + p ) * np.log( 1.0 + nu**(-1) * np.dot( np.dot( (x - mu), np.linalg.inv(Sigma) ), (x - mu) ) )

    return part1 + part2 + part3 + part4

 def _ewens_theta_likelihood (self, theta):
     '''
     returns the likelihood of theta for a given dataset
     
     '''
     if theta < 0:
         return float ('-inf')
     return self.community.S * log(theta) + lgamma(theta) - lgamma(theta + self.community.J)

def p(n,m,p):
    """Probability of m success out of n events, where an individual event succeeds with probability P... Useful for calculating <H^hat(j|w)> in semantic_information fun\
ction below"""
    try:
        return exp(lgamma(n+1)  - lgamma(n-m+1) - lgamma(m+1) +  m*log(p) + (n-m)*log(1.0-p))
    except:
        print "WARNING: domain range errer...returning 0"
        return 0

def log_Beta(alphas):
    """
    the beta function is a product-of-gammas over a gamma-of-sum
    the gamma function generalizes factorials
    tested against wolfram alpha
    """
    #return product(map(gamma,alphas)) / gamma(sum(alphas))
    return sum(lgamma(alpha) for alpha in alphas) - lgamma(sum(alphas))

def Bernstein(n, k):
    """Bernstein polynomial.

    """
    # binom
    coeff = exp(lgamma(1+n)-lgamma(1+k)-lgamma(1+n-k))

    return lambda x: coeff*x**k*(1-x)**(n-k)