Python scipy.special.comb() Examples

def get_celeba_task_pool(self, attributes, order=3, print_partition=None):
        """
        Produces partitions: a list of dictionaries (key: 0 or 1, value: list of data indices), which is
        compatible with the other methods of this class.
        """
        num_pools = 0
        partitions = []
        from scipy.special import comb
        for attr_comb in tqdm(combinations(range(attributes.shape[1]), order), desc='get_task_pool', total=comb(attributes.shape[1], order)):
            for booleans in product(range(2), repeat=order-1):
                booleans = (0,) + booleans  # only the half of the cartesian products that start with 0
                positive = np.where(np.all([attributes[:, attr] == i_booleans for (attr, i_booleans) in zip(attr_comb, booleans)], axis=0))[0]
                if len(positive) < self.num_samples_per_class:
                    continue
                negative = np.where(np.all([attributes[:, attr] == 1 - i_booleans for (attr, i_booleans) in zip(attr_comb, booleans)], axis=0))[0]
                if len(negative) < self.num_samples_per_class:
                    continue
                # inner_pool[booleans] = {0: list(negative), 1: list(positive)}
                partitions.append({0: list(negative), 1: list(positive)})
                num_pools += 1
                if num_pools == print_partition:
                    print(attr_comb, booleans)
        print('Generated {} task pools by using {} attributes from {} per pool'.format(num_pools, order, attributes.shape[1]))
        return partitions 

def NumRegressors(npix, pld_order, cross_terms=True):
    '''
    Return the number of regressors for `npix` pixels
    and PLD order `pld_order`.

    :param bool cross_terms: Include pixel cross-terms? Default :py:obj:`True`

    '''

    res = 0
    for k in range(1, pld_order + 1):
        if cross_terms:
            res += comb(npix + k - 1, k)
        else:
            res += npix
    return int(res) 

def _exact_p_value(self):
        r"""
        Computes the exact p-value of the McNemar test.

        Returns
        -------
        p_value : float
            The calculated exact p-value.

        Notes
        -----
        The one-sided exact p-value is defined as the following:

        .. math::

            p_{exact} = \sum^n_{i=n_{12}} \binom{n}{i} \frac{1}{2}^i (1 - \frac{1}{2})^{n - i}

        """
        i = self.table[0, 1]
        n = self.table[1, 0] + self.table[0, 1]
        i_n = np.arange(i + 1, n + 1)

        p_value = 1 - np.sum(comb(n, i_n) * 0.5 ** i_n * (1 - 0.5) ** (n - i_n))

        return p_value * 2 

def _construct_collection(
        orders,
        dist,
        x_lookup,
        w_lookup,
):
    """Create a collection of {abscissa: weight} key-value pairs."""
    order = numpy.min(orders)
    skew = orders-order

    # Indices and coefficients used in the calculations
    indices = numpoly.glexindex(
        order-len(dist)+1, order+1, dimensions=len(dist))
    coeffs = numpy.sum(indices, -1)
    coeffs = (2*((order-coeffs+1) % 2)-1)*comb(len(dist)-1, order-coeffs)

    collection = defaultdict(float)
    for bidx, coeff in zip(indices+skew, coeffs.tolist()):
        abscissas = [value[idx] for idx, value in zip(bidx, x_lookup)]
        weights = [value[idx] for idx, value in zip(bidx, w_lookup)]
        for abscissa, weight in zip(product(*abscissas), product(*weights)):
            collection[abscissa] += numpy.prod(weight)*coeff

    return collection 

def get_n_boards_LUT(self):
        _c = self.get_n_cards_dealt_in_transition_to_LUT()
        return {
            r: comb(N=self.rules.N_RANKS * self.rules.N_SUITS, k=_c[r], exact=True, repetition=False)
            for r in self.rules.ALL_ROUNDS_LIST
        } 

def pmf(self, input_values, x):
        """
        Calculates the probability mass function at point x.

        Parameters
        ----------
        input_values: list
            List of input parameters, in the same order as specified in the InputConnector passed to the init function
        x: list
            The point at which the pmf should be evaluated.

        Returns
        -------
        Float
            The evaluated pmf at point x.
        """

        # If the provided point is not an integer, it is converted to one
        x = int(x)
        n = input_values[0]
        p = input_values[1]
        if(x>n):
            pmf = 0
        else:
            pmf = comb(n,x)*pow(p,x)*pow((1-p),(n-x))
        self.calculated_pmf = pmf
        return pmf 

def precomputeTerms(P, ind0):
    """
    It precomputes the terms for P2M and M2M computation.

    Arguments
    ----------
    P   : int, order of the Taylor expansion.
    ind0: class, it contains the indices related to the treecode computation.
    """

    # Precompute terms for
    ind0.combII = numpy.array([], dtype=numpy.int32)
    ind0.combJJ = numpy.array([], dtype=numpy.int32)
    ind0.combKK = numpy.array([], dtype=numpy.int32)
    ind0.IImii = numpy.array([], dtype=numpy.int32)
    ind0.JJmjj = numpy.array([], dtype=numpy.int32)
    ind0.KKmkk = numpy.array([], dtype=numpy.int32)
    ind0.index_small = numpy.array([], dtype=numpy.int32)
    ind0.index_ptr = numpy.zeros(len(ind0.II) + 1, dtype=numpy.int32)
    for i in range(len(ind0.II)):
        ii, jj, kk = numpy.mgrid[0:ind0.II[i] + 1:1, 0:ind0.JJ[i] + 1:1, 0:
                                 ind0.KK[i] + 1:1].astype(numpy.int32)
        ii, jj, kk = ii.ravel(), jj.ravel(), kk.ravel()
        index_aux = numpy.zeros(len(ii), numpy.int32)
        getIndex_arr(P, len(ii), index_aux, ii, jj, kk)
        ind0.index_small = numpy.append(ind0.index_small, index_aux)
        ind0.index_ptr[i + 1] = len(index_aux) + ind0.index_ptr[i]
        ind0.combII = numpy.append(ind0.combII, comb(ind0.II[i], ii))
        ind0.combJJ = numpy.append(ind0.combJJ, comb(ind0.JJ[i], jj))
        ind0.combKK = numpy.append(ind0.combKK, comb(ind0.KK[i], kk))
        ind0.IImii = numpy.append(ind0.IImii, ind0.II[i] - ii)
        ind0.JJmjj = numpy.append(ind0.JJmjj, ind0.JJ[i] - jj)
        ind0.KKmkk = numpy.append(ind0.KKmkk, ind0.KK[i] - kk) 

def rand_index_score(clusters, classes):
    tp_plus_fp = comb(np.bincount(clusters), 2).sum()
    tp_plus_fn = comb(np.bincount(classes), 2).sum()
    A = np.c_[(clusters, classes)]
    tp = sum(comb(np.bincount(A[A[:, 0] == i, 1]), 2).sum()
             for i in set(clusters))
    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 from_bernstein_basis(cls, bp, extrapolate=None):
        """
        Construct a piecewise polynomial in the power basis
        from a polynomial in Bernstein basis.

        Parameters
        ----------
        bp : BPoly
            A Bernstein basis polynomial, as created by BPoly
        extrapolate : bool or 'periodic', optional
            If bool, determines whether to extrapolate to out-of-bounds points
            based on first and last intervals, or to return NaNs.
            If 'periodic', periodic extrapolation is used. Default is True.
        """
        dx = np.diff(bp.x)
        k = bp.c.shape[0] - 1  # polynomial order

        rest = (None,)*(bp.c.ndim-2)

        c = np.zeros_like(bp.c)
        for a in range(k+1):
            factor = (-1)**a * comb(k, a) * bp.c[a]
            for s in range(a, k+1):
                val = comb(k-a, s-a) * (-1)**s
                c[k-s] += factor * val / dx[(slice(None),)+rest]**s

        if extrapolate is None:
            extrapolate = bp.extrapolate

        return cls.construct_fast(c, bp.x, extrapolate, bp.axis) 

def from_power_basis(cls, pp, extrapolate=None):
        """
        Construct a piecewise polynomial in Bernstein basis
        from a power basis polynomial.

        Parameters
        ----------
        pp : PPoly
            A piecewise polynomial in the power basis
        extrapolate : bool or 'periodic', optional
            If bool, determines whether to extrapolate to out-of-bounds points
            based on first and last intervals, or to return NaNs.
            If 'periodic', periodic extrapolation is used. Default is True.
        """
        dx = np.diff(pp.x)
        k = pp.c.shape[0] - 1   # polynomial order

        rest = (None,)*(pp.c.ndim-2)

        c = np.zeros_like(pp.c)
        for a in range(k+1):
            factor = pp.c[a] / comb(k, k-a) * dx[(slice(None),)+rest]**(k-a)
            for j in range(k-a, k+1):
                c[j] += factor * comb(j, k-a)

        if extrapolate is None:
            extrapolate = pp.extrapolate

        return cls.construct_fast(c, pp.x, extrapolate, pp.axis) 

def _munp(self, n, c):
        def __munp(n, c):
            val = 0.0
            k = np.arange(0, n + 1)
            for ki, cnk in zip(k, sc.comb(n, k)):
                val = val + cnk * (-1) ** ki / (1.0 - c * ki)
            return np.where(c * n < 1, val * (-1.0 / c) ** n, np.inf)
        return _lazywhere(c != 0, (c,),
                          lambda c: __munp(n, c),
                          sc.gamma(n + 1)) 

def _munp(self, n, c):
        def __munp(n, c):
            val = 0.0
            k = np.arange(0, n + 1)
            for ki, cnk in zip(k, sc.comb(n, k)):
                val = val + cnk * (-1) ** ki / (1.0 - c * ki)
            return np.where(c * n < 1, val * (-1.0 / c) ** n, np.inf)
        return _lazywhere(c != 0, (c,),
                          lambda c: __munp(n, c),
                          sc.gamma(n + 1)) 

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

def _inverse_phi(self, u_loc, theta, order):
        @lru_cache(None)
        def iphi(n):
            if n:
                return sum(special.comb(n, i-1)*iphi(n-i)*sigma(i)
                           for i in range(1, n+1))
            return numpy.e**(-u_loc**(1/theta))
        @lru_cache(None)
        def sigma(n):
            return self._sigma(u_loc, theta, n)
        return iphi(order) 

def _mom(self, k, C, Ci, loc):
        scale = numpy.dot(C, C.T)
        out = 0.
        for idx, kdx in enumerate(numpy.ndindex(*[_+1 for _ in k])):
            coef = numpy.prod(special.comb(k.T, kdx).T, 0)
            diff = k.T - kdx
            pos = diff >= 0
            diff = diff*pos
            pos = numpy.all(pos)
            loc_ = numpy.prod(loc**diff)

            out += pos*coef*loc_*isserlis_moment(tuple(kdx), scale)

        return float(out) 

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 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 ser2ber(q,n,d,t,ps):
    """
    Converts symbol error rate to bit error rate. Taken from Ziemer and
    Tranter page 650. Necessary when comparing different types of block codes.
    
    parameters
    ----------  
    q: size of the code alphabet for given modulation type (BPSK=2)
    n: number of channel bits
    d: distance (2e+1) where e is the number of correctable errors per code word.
       For hamming codes, e=1, so d=3.
    t: number of correctable errors per code word
    ps: symbol error probability vector
    
    returns
    -------
    ber: bit error rate
    
    """
    lnps = len(ps) # len of error vector
    ber = np.zeros(lnps) # inialize output vector
    for k in range(0,lnps): # iterate error vector
        ser = ps[k] # channel symbol error rate
        sum1 = 0 # initialize sums
        sum2 = 0
        for i in range(t+1,d+1):
            term = special.comb(n,i)*(ser**i)*((1-ser))**(n-i)
            sum1 = sum1 + term
        for i in range(d+1,n+1):
            term = (i)*special.comb(n,i)*(ser**i)*((1-ser)**(n-i))
            sum2 = sum2+term
        ber[k] = (q/(2*(q-1)))*((d/n)*sum1+(1/n)*sum2)
    
    return ber 

def my_althom_likelihod(a,b):
	baseq_major=[float(i) for i in a.split(',')[:-1]]
	baseq_minor=[float(i) for i in b.split(',')[:-1]]
	depth=len(baseq_major)+len(baseq_minor)
	alt=len(baseq_minor)
	q=math.log10(1)
	q=sum(math.log10(1-0.1**(i/10)) for i in baseq_minor)
	q=q+sum(math.log10(0.1**(i/10)) for i in baseq_major)
	l=math.log10(comb(depth,alt,exact=True))+q
	return(10**l) 

def my_refhom_likelihod(a,b):
	baseq_major=[float(i) for i in a.split(',')[:-1]]
	baseq_minor=[float(i) for i in b.split(',')[:-1]]
	depth=len(baseq_major)+len(baseq_minor)
	alt=len(baseq_minor)
	q=math.log10(1)
	q=sum(math.log10(1-0.1**(i/10)) for i in baseq_major)
	q=q+sum(math.log10(0.1**(i/10)) for i in baseq_minor)
	l=math.log10(comb(depth,alt,exact=True))+q
	#return(10**q)	
	return(10**l) 

def my_het_likelihood(a,b,c,d):
	depth=sum([int(a),int(b),int(c),int(d)])
	alt=sum([int(c),int(d)])
	l=math.log10(comb(depth,alt,exact=True))+math.log10(0.5)*depth
	l=10**l
	return(l)
#	math.log10(comb(2000,1000,exact=True))+math.log10(0.5)*2000
#	return(0.5**depth) 

def my_mosaic_likelihood(a,b,c,d,e,f):
	depth=sum([int(a),int(b),int(c),int(d)])
	alt=sum([int(c),int(d)])
	r=0
	baseq_major=[float(i) for i in e.split(',')[:-1]]
	baseq_minor=[float(i) for i in f.split(',')[:-1]]
	r=sum([0.1**(float(i)/10) for i in baseq_major])
	r=r+sum([1-0.1**(float(i)/10) for i in baseq_minor])
	l=beta(r+1, depth-r+1)
	if l >0:
		l=math.log10(l)+math.log10(comb(depth,alt,exact=True))
		l=10**l
##	return(float(beta(r+1, depth-r+1)))
	return(l) 

def __len__(self):
        if self.length is not None:
            return self.length
        return int(len(self.keys) * (comb(self.max_samples, self.batch_k) + comb(self.min_samples, self.batch_k))/2) 

def confusion_index(v1, v2):
    cmatrix = contingency(v1, v2)
    size = np.size(v1)
    sum_rows = np.sum(cmatrix, 0)
    sum_cols = np.sum(cmatrix, 1)
    N = comb(size, 2)
    TP = np.sum(list(map(lambda x: comb(x, 2), cmatrix)))
    FN = np.sum(list(map(lambda x: comb(x, 2), sum_rows))) - TP
    FP = np.sum(list(map(lambda x: comb(x, 2), sum_cols))) - TP
    TN = N - TP - FN - FP
    return TP, FN, FP, TN 

def test_prior():
    K = 10
    T = 100

    es = EventSegment(K)
    mp = es.model_prior(T)[0]

    p_bound = np.zeros((T, K-1))
    norm = comb(T-1, K-1)
    for t in range(T-1):
        for k in range(K-1):
            # See supplementary material of Neuron paper
            # https://doi.org/10.1016/j.neuron.2017.06.041
            p_bound[t+1, k] = comb(t, k) * comb(T-t-2, K-k-2) / norm
    p_bound = np.cumsum(p_bound, axis=0)

    mp_gt = np.zeros((T, K))
    for k in range(K):
        if k == 0:
            mp_gt[:, k] = 1 - p_bound[:, 0]
        elif k == K - 1:
            mp_gt[:, k] = p_bound[:, k-1]
        else:
            mp_gt[:, k] = p_bound[:, k-1] - p_bound[:, k]

    assert np.all(np.isclose(mp, mp_gt)),\
        "Prior does not match analytic solution" 

def test_comb_zeros(self):
        assert_equal(special.comb(2, 3, exact=True), 0)
        assert_equal(special.comb(-1, 3, exact=True), 0)
        assert_equal(special.comb(2, -1, exact=True), 0)
        assert_equal(special.comb(2, -1, exact=False), 0)
        assert_array_almost_equal(special.comb([2, -1, 2, 10], [3, 3, -1, 3]),
                [0., 0., 0., 120.]) 

def test_comb_with_np_int64(self):
        n = 70
        k = 30
        np_n = np.int64(n)
        np_k = np.int64(k)
        assert_equal(special.comb(np_n, np_k, exact=True),
                     special.comb(n, k, exact=True)) 

def test_comb(self):
        assert_array_almost_equal(special.comb([10, 10], [3, 4]), [120., 210.])
        assert_almost_equal(special.comb(10, 3), 120.)
        assert_equal(special.comb(10, 3, exact=True), 120)
        assert_equal(special.comb(10, 3, exact=True, repetition=True), 220)

        assert_allclose([special.comb(20, k, exact=True) for k in range(21)],
                        special.comb(20, list(range(21))), atol=1e-15)

        ii = np.iinfo(int).max + 1
        assert_equal(special.comb(ii, ii-1, exact=True), ii)

        expected = 100891344545564193334812497256
        assert_equal(special.comb(100, 50, exact=True), expected) 

def _nCr(n, r):
    """Number of combinations of r items out of a set of n.  Equals n!/(r!(n-r)!)"""
    #f = _math.factorial; return f(n) / f(r) / f(n-r)
    return _spspecial.comb(n, r) 

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