Python scipy.misc.comb() Examples

def mnc2mc(mnc, wmean = True):
    '''convert non-central to central moments, uses recursive formula
    optionally adjusts first moment to return mean

    '''
    n = len(mnc)
    mean = mnc[0]
    mnc = [1] + list(mnc)    # add zero moment = 1
    mu = [] #np.zeros(n+1)
    for n,m in enumerate(mnc):
        mu.append(0)
        #[comb(n-1,k,exact=1) for k in range(n)]
        for k in range(n+1):
            mu[n] += (-1)**(n-k) * comb(n,k,exact=1) * mnc[k] * mean**(n-k)
    if wmean:
        mu[1] = mean
    return mu[1:] 

def index(self, text, terms=None, **kwargs):

        """
        Index all term pair distances.

        Args:
            text (Text): The source text.
            terms (list): Terms to index.
        """

        self.clear()

        # By default, use all terms.
        terms = terms or text.terms.keys()

        pairs = combinations(terms, 2)
        count = comb(len(terms), 2)

        for t1, t2 in bar(pairs, expected_size=count, every=1000):

            # Set the Bray-Curtis distance.
            score = text.score_braycurtis(t1, t2, **kwargs)
            self.set_pair(t1, t2, score) 

def mc2mnc(mc):
    '''convert central to non-central moments, uses recursive formula
    optionally adjusts first moment to return mean

    '''
    n = len(mc)
    mean = mc[0]
    mc = [1] + list(mc)    # add zero moment = 1
    mc[1] = 0  # define central mean as zero for formula
    mnc = [1, mean] # zero and first raw moments
    for nn,m in enumerate(mc[2:]):
        n=nn+2
        mnc.append(0)
        for k in range(n+1):
            mnc[n] += comb(n,k,exact=1) * mc[k] * mean**(n-k)

    return mnc[1:] 

def mnc2mc(mnc, wmean = True):
    '''convert non-central to central moments, uses recursive formula
    optionally adjusts first moment to return mean

    '''
    n = len(mnc)
    mean = mnc[0]
    mnc = [1] + list(mnc)    # add zero moment = 1
    mu = [] #np.zeros(n+1)
    for n,m in enumerate(mnc):
        mu.append(0)
        #[comb(n-1,k,exact=1) for k in range(n)]
        for k in range(n+1):
            mu[n] += (-1)**(n-k) * comb(n,k,exact=1) * mnc[k] * mean**(n-k)
    if wmean:
        mu[1] = mean
    return mu[1:] 

def cum2mc(kappa):
    '''convert non-central moments to cumulants
    recursive formula produces as many cumulants as moments

    References
    ----------
    Kenneth Lange: Numerical Analysis for Statisticians, page 40
    (http://books.google.ca/books?id=gm7kwttyRT0C&pg=PA40&lpg=PA40&dq=convert+cumulants+to+moments&source=web&ots=qyIaY6oaWH&sig=cShTDWl-YrWAzV7NlcMTRQV6y0A&hl=en&sa=X&oi=book_result&resnum=1&ct=result)


    '''
    mc = [1,0.0] #_kappa[0]]  #insert 0-moment and mean
    kappa0 = kappa[0]
    kappa = [1] + list(kappa)
    for nn,m in enumerate(kappa[2:]):
        n = nn+2
        mc.append(0)
        for k in range(n-1):
            mc[n] += comb(n-1,k,exact=1) * kappa[n-k]*mc[k]

    mc[1] = kappa0 # insert mean as first moments by convention
    return mc[1:] 

def kendalls_tau(X):
    """
    Calculates a generalized Kendall's tau for a data set given by X, as 
    described by "Multivariate Extensions of Spearman's Rho and Related Statistics"
    
    Inputs:
      X - the input data, should be a numpy array of shape = M x N, where
          M is the number of samples, and N is the dimensionality of the data
    """
    M = X.shape[0]
    N = X.shape[1]
    if N<2:
        raise ValueError('To calculate Kendall\'s Tau, need data of dimensionality >= 2')
    
    ktau = 0.0
    for dim1 in range(0,N-1):
        for dim2 in range(dim1+1,N):
            (t,p) = kendalltau(X[:,dim1],X[:,dim2])
            ktau = ktau + t
    # normalize
    ktau = ktau / comb(N,2)
    return ktau 

def spearmans_rho(X):
    """
    Calculates a generalized Spearman's rho for a data set given by X, as 
    described by "Multivariate Extensions of Spearman's Rho and Related Statistics"
    Inputs:
      X - the input data, should be a numpy array of shape = M x N, where
          M is the number of samples, and N is the dimensionality of the data
    """
    M = X.shape[0]
    N = X.shape[1]
    if N<2:
        raise ValueError('To calculate Spearman\'s Rho, need data of dimensionality >= 2')
    
    srho = 0.0
    for dim1 in range(0,N-1):
        for dim2 in range(dim1+1,N):
            (r,p) = spearmanr(X[:,dim1],X[:,dim2])
            srho = srho + r
    # normalize
    srho = srho / comb(N,2)
    return srho 

def cal_pvalue(c00, c11, c01, c10):
    """
    c00: int, argeed number on positive
    c01, c10: int, disagreed number
    c11: int, agreed number on negative
    return: p value
    """

    b = min(c01, c10)
    n = c01 + c10
    p = 0
    for i in range(b, n+1):
        cur_p = 1
        for j in range(i):
            cur_p *= ((n - j) / (2.0 * (i - j)))
        if not np.isinf(cur_p):# and n-i < 1000: 
            # print(n-i, cur_p)
            for i in range(n-i):
                cur_p /= 2.0
            p += cur_p
            # print(cur_p, p)
        # p = p + comb(c01+c10, j) / denominator
    # p /= (2 ** (c01 + c10 - 1))
    return 2 * p 

def generate_hamming_distance_payoff_distribution(max_sent_len, vocab_size, tau=1.):
    """compute the q distribution for Hamming Distance (substitution only) as in the RAML paper"""
    probs = dict()
    Z_qs = dict()
    for sent_len in range(1, max_sent_len + 1):
        counts = [1.]  # e = 0, count = 1
        for e in range(1, sent_len + 1):
            # apply the rescaling trick as in https://gist.github.com/norouzi/8c4d244922fa052fa8ec18d8af52d366
            count = comb(sent_len, e) * math.exp(-e / tau) * ((vocab_size - 1) ** (e - e / tau))
            counts.append(count)

        Z_qs[sent_len] = Z_q = sum(counts)
        prob = [count / Z_q for count in counts]
        probs[sent_len] = prob

        # print('sent_len=%d, %s' % (sent_len, prob))

    return probs, Z_qs 

def ab_filter_exp_cascade(self,tau, n, step = 0.001):
    """ create an ExpCascade filter and return arrays for the coefficients """
    from scipy.misc import comb
    tau = float(tau)
    n = int(n)
    if tau < 0:
        raise Exception("Negative time constants are not implemented")
    if n < 0:
        raise Exception("Negative cascade number is not allowed")
    if tau == 0:
        a = np.array([1.0])
        b = np.array([])
        return [a,b]
    tauC = tau/float(n) if n > 0 else tau
    c = np.exp(-step/tauC)
    N = n + 1
    a = np.array([(-c)**(i)*comb(N,i) for i in range(N+1)])
    b = np.array([(1.0-c)**N])
    return [a,b] 

def mc2mnc(mc):
    '''convert central to non-central moments, uses recursive formula
    optionally adjusts first moment to return mean

    '''
    n = len(mc)
    mean = mc[0]
    mc = [1] + list(mc)    # add zero moment = 1
    mc[1] = 0  # define central mean as zero for formula
    mnc = [1, mean] # zero and first raw moments
    for nn,m in enumerate(mc[2:]):
        n=nn+2
        mnc.append(0)
        for k in range(n+1):
            mnc[n] += comb(n,k,exact=1) * mc[k] * mean**(n-k)

    return mnc[1:] 

def cum2mc(kappa):
    '''convert non-central moments to cumulants
    recursive formula produces as many cumulants as moments

    References
    ----------
    Kenneth Lange: Numerical Analysis for Statisticians, page 40
    (http://books.google.ca/books?id=gm7kwttyRT0C&pg=PA40&lpg=PA40&dq=convert+cumulants+to+moments&source=web&ots=qyIaY6oaWH&sig=cShTDWl-YrWAzV7NlcMTRQV6y0A&hl=en&sa=X&oi=book_result&resnum=1&ct=result)


    '''
    mc = [1,0.0] #_kappa[0]]  #insert 0-moment and mean
    kappa0 = kappa[0]
    kappa = [1] + list(kappa)
    for nn,m in enumerate(kappa[2:]):
        n = nn+2
        mc.append(0)
        for k in range(n-1):
            mc[n] += comb(n-1,k,exact=1) * kappa[n-k]*mc[k]

    mc[1] = kappa0 # insert mean as first moments by convention
    return mc[1:] 

def _munp(self, n, c):
        k = np.arange(0,n+1)
        val = (1.0/c)**n * np.sum(comb(n,k)*(-1)**k / (1.0+c*k),axis=0)
        return where(c*n > -1, val, inf) 

def pdf(self, x, k, n, p):
        '''distribution of success runs of length k or more

        Parameters
        ----------
        x : float
            count of runs of length n
        k : int
            length of runs
        n : int
            total number of observations or trials
        p : float
            probability of success in each Bernoulli trial

        Returns
        -------
        pdf : float
            probability that x runs of length of k are observed

        Notes
        -----
        not yet vectorized

        References
        ----------
        Muselli 1996, theorem 3
        '''

        q = 1-p
        m = np.arange(x, (n+1)//(k+1)+1)[:,None]
        terms = (-1)**(m-x) * comb(m, x) * p**(m*k) * q**(m-1) \
                * (comb(n - m*k, m - 1) + q * comb(n - m*k, m))
        return terms.sum(0) 

def runs_prob_odd(self, r):
        n0, n1 = self.n0, self.n1
        k = (r+1)//2
        tmp0 = comb(n0-1, k-1)
        tmp1 = comb(n1-1, k-2)
        tmp3 = comb(n0-1, k-2)
        tmp4 = comb(n1-1, k-1)
        return (tmp0 * tmp1 + tmp3 * tmp4)  / self.comball 

def __init__(self, n0, n1):
        self.n0 = n0
        self.n1 = n1
        self.n = n = n0 + n1
        self.comball = comb(n, n1) 

def mnc2cum(mnc):
    '''convert non-central moments to cumulants
    recursive formula produces as many cumulants as moments

    http://en.wikipedia.org/wiki/Cumulant#Cumulants_and_moments
    '''
    mnc = [1] + list(mnc)
    kappa = [1]
    for nn,m in enumerate(mnc[1:]):
        n = nn+1
        kappa.append(m)
        for k in range(1,n):
            kappa[n] -= comb(n-1,k-1,exact=1) * kappa[k]*mnc[n-k]

    return kappa[1:] 

def resolve_conflicts(meta_all, snames):
    fcutoff = 50
    # meta_all['nconflict'] = meta_all['conflict_with'].apply(lambda l: np.sum(l))
    # meta_all = meta_all.sort_values('nconflict', ascending=False)
    conf_pairs = meta_all[meta_all['conflict_with'].apply(bool)]['conflict_with']
    conf_list = [[list(conf_pairs.index).index(c) for c in p] for p in conf_pairs]
    nconf = len(conf_list)
    conf_tree = nbtree.NBTree(nconf)
    remv_exp_last = 0
    print ("processing conflict list with length {}".format(nconf))
    for level in range(fcutoff, nconf):
        t = time.time()
        comb_exp = np.sum([misc.comb(level, posslen) for posslen in range(fcutoff, level + 1)])
        remv_exp = comb_exp - remv_exp_last*2
        remv_exp_last = comb_exp
        nremoved = 0
        nsearched = 0
        print ("processing level: " + str(level))
        print ("expected removal: " + str(remv_exp))
        for nd in conf_tree.get_nodes(level=level, unpack=True):
            print ("searched: {0}, removed: {1}".format(nsearched, nremoved), end='\r')
            if nremoved >= remv_exp:
                break
            if nd % 2:
                nsearched += 1
                fal = np.array(conf_tree.path_to_node(nd)) % 2
                if np.sum(fal) >= fcutoff:
                    conf_tree.remove_subtree(nd)
                    nremoved += 1
        print ("process time: {} s".format(time.time() - t))
    for conid1, pair in enumerate(conf_list):
        for conid2 in pair:
            conidl = min(conid1, conid2)
            conidh = max(conid1, conid2)
            print ("processing pair:" + str((conidl, conidh)))
            for ndl in conf_tree.get_nodes(level=conidl + 1, unpack=True):
                if not ndl / 2:
                    for ndh in conf_tree.children(ndl).get_nodes(level=conidh + 1, unpack=True):
                        if not ndh / 2:
                            conf_tree.remove_subtree(ndh)
    return conf_tree 

def __init__(self, name, layout_constraint, entity_constraint, 
                             orig_layout_constraint=None, orig_entity_constraint=None, 
                             sample_new_num_count=None, is_pg=False):
        super(Layout, self).__init__(name, level="Layout", node_type="and", is_pg=is_pg)
        self.layout_constraint = layout_constraint
        self.entity_constraint = entity_constraint
        self.number = Number(min_level=layout_constraint["Number"][0], max_level=layout_constraint["Number"][1])
        self.position = Position(pos_type=layout_constraint["Position"][0], pos_list=layout_constraint["Position"][1])
        self.uniformity = Uniformity(min_level=layout_constraint["Uni"][0], max_level=layout_constraint["Uni"][1])
        self.number.sample()
        self.position.sample(self.number.get_value())
        self.uniformity.sample()
        # store initial layout_constraint and entity_constraint for answer generation
        if orig_layout_constraint is None:
            self.orig_layout_constraint = copy.deepcopy(self.layout_constraint)
        else:
            self.orig_layout_constraint = orig_layout_constraint
        if orig_entity_constraint is None:
            self.orig_entity_constraint = copy.deepcopy(self.entity_constraint)
        else:
            self.orig_entity_constraint = orig_entity_constraint
        if sample_new_num_count is None:
            self.sample_new_num_count = dict()
            most_num = len(self.position.values)
            for i in range(layout_constraint["Number"][0], layout_constraint["Number"][1] + 1):
                self.sample_new_num_count[i] = [comb(most_num, i + 1), []]
        else:
            self.sample_new_num_count = sample_new_num_count 

def bernstein_poly(i, n, t):
    """
     The Bernstein polynomial of n, i as a function of t
    """

    return comb(n, i) * ( t**(n-i) ) * (1 - t)**i 

def rand_index_score(label_pred, label_gt):
    tp_plus_fp = comb(np.bincount(label_pred), 2).sum()
    tp_plus_fn = comb(np.bincount(label_gt), 2).sum()
    A = np.c_[(label_pred, label_gt)]
    tp = sum(comb(np.bincount(A[A[:, 0] == i, 1]), 2).sum()
             for i in set(label_pred))
    fp = tp_plus_fp - tp
    fn = tp_plus_fn - tp
    tn = comb(len(A), 2) - tp - fp - fn
    return (tp + tn) / (tp + fp + fn + tn) 

def graphlets_counts_4(self, n, k):
        res = np.array([comb(k,4), # cliques
                        comb(k,3)*(n-k), #triangle + dot
                        comb(k,2)*comb(n-k,2), #two dots and edge
                        comb(n-k,4) + comb(n-k,3)*k]) # empty
        return (res/comb(n,4)) 

def graphlets_counts_3(self, n,k):
        res = np.array([comb(n-k,3) + comb(n-k,2)*k, #empty
                        comb(k,2)*(n-k), #one edge
                        0.0, # two edges
                        comb(k,3)]) #triangles
        return (res/comb(n,3)) 

def bernstein_poly(i, n, t):
    """
     The Bernstein polynomial of n, i as a function of t
    """

    return comb(n, i) * ( t**(n-i) ) * (1 - t)**i 

def test_big(self):
        p = pascal(50)
        assert_equal(p[-1, -1], comb(98, 49, exact=True)) 

def _kolmog1(x,n):
    if x <= 0:
        return 0
    if x >= 1:
        return 1
    j = np.arange(np.floor(n*(1-x))+1)
    return 1 - x * np.sum(np.exp(np.log(misc.comb(n,j))
                                       + (n-j) * np.log(1-x-j/float(n))
                                       + (j-1) * np.log(x+j/float(n)))) 

def bilinear(b, a, fs=1.0):
    """Return a digital filter from an analog one using a bilinear transform.

    The bilinear transform substitutes ``(z-1) / (z+1)`` for ``s``.
    """
    fs = float(fs)
    a, b = map(atleast_1d, (a, b))
    D = len(a) - 1
    N = len(b) - 1
    artype = float
    M = max([N, D])
    Np = M
    Dp = M
    bprime = numpy.zeros(Np + 1, artype)
    aprime = numpy.zeros(Dp + 1, artype)
    for j in range(Np + 1):
        val = 0.0
        for i in range(N + 1):
            for k in range(i + 1):
                for l in range(M - i + 1):
                    if k + l == j:
                        val += (comb(i, k) * comb(M - i, l) * b[N - i] *
                                pow(2 * fs, i) * (-1) ** k)
        bprime[j] = real(val)
    for j in range(Dp + 1):
        val = 0.0
        for i in range(D + 1):
            for k in range(i + 1):
                for l in range(M - i + 1):
                    if k + l == j:
                        val += (comb(i, k) * comb(M - i, l) * a[D - i] *
                                pow(2 * fs, i) * (-1) ** k)
        aprime[j] = real(val)

    return normalize(bprime, aprime) 

def lp2bs(b, a, wo=1.0, bw=1.0):
    """
    Transform a lowpass filter prototype to a highpass filter.

    Return an analog band-stop filter with center frequency `wo` and
    bandwidth `bw` from an analog low-pass filter prototype with unity
    cutoff frequency, in transfer function ('ba') representation.

    """
    a, b = map(atleast_1d, (a, b))
    D = len(a) - 1
    N = len(b) - 1
    artype = mintypecode((a, b))
    M = max([N, D])
    Np = M + M
    Dp = M + M
    bprime = numpy.zeros(Np + 1, artype)
    aprime = numpy.zeros(Dp + 1, artype)
    wosq = wo * wo
    for j in range(Np + 1):
        val = 0.0
        for i in range(0, N + 1):
            for k in range(0, M - i + 1):
                if i + 2 * k == j:
                    val += (comb(M - i, k) * b[N - i] *
                            (wosq) ** (M - i - k) * bw ** i)
        bprime[Np - j] = val
    for j in range(Dp + 1):
        val = 0.0
        for i in range(0, D + 1):
            for k in range(0, M - i + 1):
                if i + 2 * k == j:
                    val += (comb(M - i, k) * a[D - i] *
                            (wosq) ** (M - i - k) * bw ** i)
        aprime[Dp - j] = val

    return normalize(bprime, aprime) 

def lp2bp(b, a, wo=1.0, bw=1.0):
    """
    Transform a lowpass filter prototype to a bandpass filter.

    Return an analog band-pass filter with center frequency `wo` and
    bandwidth `bw` from an analog low-pass filter prototype with unity
    cutoff frequency, in transfer function ('ba') representation.

    """
    a, b = map(atleast_1d, (a, b))
    D = len(a) - 1
    N = len(b) - 1
    artype = mintypecode((a, b))
    ma = max([N, D])
    Np = N + ma
    Dp = D + ma
    bprime = numpy.zeros(Np + 1, artype)
    aprime = numpy.zeros(Dp + 1, artype)
    wosq = wo * wo
    for j in range(Np + 1):
        val = 0.0
        for i in range(0, N + 1):
            for k in range(0, i + 1):
                if ma - i + 2 * k == j:
                    val += comb(i, k) * b[N - i] * (wosq) ** (i - k) / bw ** i
        bprime[Np - j] = val
    for j in range(Dp + 1):
        val = 0.0
        for i in range(0, D + 1):
            for k in range(0, i + 1):
                if ma - i + 2 * k == j:
                    val += comb(i, k) * a[D - i] * (wosq) ** (i - k) / bw ** i
        aprime[Dp - j] = val

    return normalize(bprime, aprime) 

def _munp(self, n, c):
        k = np.arange(0,n+1)
        val = (1.0/c)**n * np.sum(comb(n,k)*(-1)**k / (1.0+c*k),axis=0)
        return where(c*n > -1, val, inf)