def test_cross_validator_with_default_indices():
    n_samples = 4
    n_unique_labels = 4
    n_folds = 2
    p = 2
    n_iter = 10  # (the default value)

    X = np.array([[1, 2], [3, 4], [5, 6], [7, 8]])
    X_1d = np.array([1, 2, 3, 4])
    y = np.array([1, 1, 2, 2])
    labels = np.array([1, 2, 3, 4])
    loo = LeaveOneOut()
    lpo = LeavePOut(p)
    kf = KFold(n_folds)
    skf = StratifiedKFold(n_folds)
    lolo = LeaveOneLabelOut()
    lopo = LeavePLabelOut(p)
    ss = ShuffleSplit(random_state=0)
    ps = PredefinedSplit([1, 1, 2, 2])  # n_splits = np of unique folds = 2

    n_splits = [n_samples, comb(n_samples, p), n_folds, n_folds,
                n_unique_labels, comb(n_unique_labels, p), n_iter, 2]

    for i, cv in enumerate([loo, lpo, kf, skf, lolo, lopo, ss, ps]):
        # Test if get_n_splits works correctly
        assert_equal(n_splits[i], cv.get_n_splits(X, y, labels))

        # Test if the cross-validator works as expected even if
        # the data is 1d
        np.testing.assert_equal(list(cv.split(X, y, labels)),
                                list(cv.split(X_1d, y, labels)))
        # Test that train, test indices returned are integers
        for train, test in cv.split(X, y, labels):
            assert_equal(np.asarray(train).dtype.kind, 'i')
            assert_equal(np.asarray(train).dtype.kind, 'i')

def findF(k,p,q,x, side):
    
    #print k, p, q, x
    
    if k in Fdict:
        if p in Fdict[k]:
            if q in Fdict[k][p]:
                if x in Fdict[k][p][q]:
                    if side in Fdict[k][p][q][x]:
                        return Fdict[k][p][q][x][side]
    
    if q>0:
        if side==LEFT: #left
            ret=0
            for l in xrange(min(p*q+1, x+1)):
                ret+=findF(k, p, 0, x-l, LEFT)*sc.comb(p*q,l)
            
        elif side==RIGHT: #right
            ret=0
            for l in xrange(min(p*q+1, x+1)):
                if l==0:
                    ret+=0 #findF(k, p, 0, x-l, LEFT)  it requires at least one edge from left to right.
                else:
                    
                    ret+=findF(k, p, 0, x-l, LEFT)*sc.comb(q,l) #external nodes can only connect the "root" of the right tree
                    
                    #ret+=findF(k, p, 0, x-l, LEFT)*sc.comb(q,1)*sc.comb((p-1)*q,l-1)   
    else:
        if k==0:
            if p==1:
                if q==0:
                    if x==0:
                        ret=1
                    else:
                        ret=0
                else:
                    print "Error 1"
            else:
                print "Errro 0"
        else:
            ret=0
            for i in xrange(k):
                tmp=0
                for j in xrange(k-i-1,x-i):
                    tmp += findF(k-i-1,k-i,i+1,j, LEFT)*findF(i, i+1, k-i, x-j, RIGHT)
                ret += tmp*sc.comb(k-1,i)
                
    
    if k not in Fdict:
        Fdict[k]={}
    if p not in Fdict[k]:
        Fdict[k][p]={}
    if q not in Fdict[k][p]:
        Fdict[k][p][q]={}
    if x not in Fdict[k][p][q]: 
        Fdict[k][p][q][x]={}
    Fdict[k][p][q][x][side]=ret
     
    
    return ret

def brive(N, replace_zeros=True):
    """
    The brive estimator

    Parameters
    ----------
    N : np.array, int
       Counts vector
    replace_zeros: bool
       Replaces zeros with uniform prior

    Returns
    -------
    pvals
    """
    N = N.astype(np.int)
    n = sum(N)
    pvals = np.zeros(len(N), dtype=np.float64)
    for i in range(len(N)):
        if N[i]==0 or N[i]==1: continue
        trials = [comb(t-1, N[i]-1) / (t * (comb(n, N[i])))
                  for t in range(N[i], n+1)]
        pvals[i] = (float(N[i]-1)) * sum(trials)

    if replace_zeros:
        m = sum(pvals)
        if 0 < m < 1 and (pvals==0).sum() > 0:
            pvals[pvals==0] = (1 - m) / (pvals==0).sum()
    return pvals

def solution3():
    '''
    >>> solution3()
    137846528820
    '''
    from scipy.misc import comb
    print comb(40, 20, exact=True)

def rare(y, size):
	notabs = ~np.isnan(y)
	t = y[notabs]
	N = np.sum(t)
	diff = N - t
	rare_calc = np.sum(1 - comb(diff, size)/comb(N, size))
	return rare_calc

    def eval_gof( self, data):
        probs = self.eval(data.levels)
        expected = probs * data.Ntrials
        # Compute various goodness of fit statistics
        LL = -2*sum((data.Ncorr*np.log(expected/(data.Ncorr+np.finfo(float).eps))
            +(data.Ntrials-data.Ncorr)*np.log((data.Ntrials-expected)/(data.Ntrials-data.Ncorr+np.finfo(float).eps))));
        X2 = sum((data.Ncorr-expected)**2./expected/(1.-probs))

        # Adding eps to avoid log(0) Nans
        self.prinsNLL = -sum( data.Ncorr*np.log(probs+np.finfo(float).eps)+(data.Ntrials-data.Ncorr)*np.log(1.-probs+np.finfo(float).eps) )

        # Treutwein/Strasburger 1999 Eq 6 (likelihood of the data)
        L_ts = 2**(sum( data.Ntrials ))
        LL_ts = 0.0
        #L_ts = 1.0
        for level in np.arange( len(data.levels) ):
            # TODO: is right to use observed data or function values?: next two lines can chg to try fitted
            thisN = data.Ntrials[level]
            thisCorr = data.Ncorr[level]
            L_ts *= misc.comb( thisN, thisCorr ) * (probs[level]**thisCorr) * (1.0 - probs[level])**(thisN-thisCorr) 
            LL_ts += np.log(misc.comb( thisN, thisCorr )) + thisCorr*np.log(probs[level]) +np.log(1.0 - probs[level])*(thisN-thisCorr) 

        #TODO: This is how Prins' clamps the lapse.  Parameterize.
        if (self.params[self.PARAM_UPPER] < 0) or (self.params[self.PARAM_UPPER] > 0.05):
            self.prinsNLL=np.inf

        return probs,LL,X2,L_ts,LL_ts,self.prinsNLL

    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 daub(p):
    """
    The coefficients for the FIR low-pass filter producing Daubechies wavelets.

    p>=1 gives the order of the zero at f=1/2.
    There are 2p filter coefficients.

    Parameters
    ----------
    p : int
        Order of the zero at f=1/2, can have values from 1 to 34.

    """
    sqrt = np.sqrt
    if p < 1:
        raise ValueError("p must be at least 1.")
    if p==1:
        c = 1/sqrt(2)
        return np.array([c,c])
    elif p==2:
        f = sqrt(2)/8
        c = sqrt(3)
        return f*np.array([1+c,3+c,3-c,1-c])
    elif p==3:
        tmp  = 12*sqrt(10)
        z1 = 1.5 + sqrt(15+tmp)/6 - 1j*(sqrt(15)+sqrt(tmp-15))/6
        z1c = np.conj(z1)
        f = sqrt(2)/8
        d0 = np.real((1-z1)*(1-z1c))
        a0 = np.real(z1*z1c)
        a1 = 2*np.real(z1)
        return f/d0*np.array([a0, 3*a0-a1, 3*a0-3*a1+1, a0-3*a1+3, 3-a1, 1])
    elif p<35:
        # construct polynomial and factor it
        if p<35:
            P = [comb(p-1+k,k,exact=1) for k in range(p)][::-1]
            yj = np.roots(P)
        else:  # try different polynomial --- needs work
            P = [comb(p-1+k,k,exact=1)/4.0**k for k in range(p)][::-1]
            yj = np.roots(P) / 4
        # for each root, compute two z roots, select the one with |z|>1
        # Build up final polynomial
        c = np.poly1d([1,1])**p
        q = np.poly1d([1])
        for k in range(p-1):
            yval = yj[k]
            part = 2*sqrt(yval*(yval-1))
            const = 1-2*yval
            z1 = const + part
            if (abs(z1)) < 1:
                z1 = const - part
            q = q * [1,-z1]

        q = c * np.real(q)
        # Normalize result
        q = q / np.sum(q) * sqrt(2)
        return q.c[::-1]
    else:
        raise ValueError("Polynomial factorization does not work "
              "well for p too large.")

def ARI(trueLab, predLab):
    """
    compute adjusted rand index, ranges in [-1, 1] with random assignment score 0 and perfect score 1
    :param trueLab: ground truth labels
    :param predLab: predicted labels
    :return: adjusted rand index
    """
    n = len(trueLab)
    trueLab = np.array(trueLab)
    predLab = np.array(predLab)

    trueCluster = dict(zip(set(trueLab), [np.where(trueLab == x)[0] for x in set(trueLab)]))
    predCluster = dict(zip(set(predLab), [np.where(predLab == x)[0] for x in set(predLab)]))

    nTrue = len(trueCluster)
    nPred = len(predCluster)

    cTable = np.zeros((nTrue, nPred))

    for i in range(nTrue):
        for j in range(nPred):
            cTable[i, j] = len(np.intersect1d(trueCluster.values()[i], predCluster.values()[j]))

    a = comb(np.sum(cTable, axis=1), 2).sum()
    b = comb(np.sum(cTable, axis=0), 2).sum()
    c = comb(n, 2)

    return (comb(cTable, 2).sum() - (a * b) / c) / (0.5 * (a + b) - (a * b) / c)

def get_motifspace_size(q, n):
    """return length of search space according to equation which is mentioned in Section 3.1 of the paper"""
    return reduce(
        lambda x, y: x + (int(sc.comb(q, y, exact=True)) * 4 ** (q - y)),
        [i for i in range(1, n + 1)],
        int(sc.comb(q, 0, exact=True)) * 4 ** (q - 0),
    )

def estimate(E, r, f, lp):
    if f <= r:
        return 0.0
    if E - f < lp:
        return 1.0
    probs = [(1 - (1.0 * scal(comb(E - f, lp)) / scal(comb(e, lp)))) for e in range(E, E - (r+ 1), -1)]
    return reduce(lambda x, y: x * y, probs, 1.0)

def bilinear(b, a, fs=1.0):
    """Return a digital filter from an analog filter using the bilinear transform.
    The bilinear transform substitutes ``(z-1) / (z+1``) for ``s``.
    """
    #This function has been copied out of scipy
    fs =float(fs)
    a,b = map(num.atleast_1d,(a,b))
    D = len(a) - 1
    N = len(b) - 1
    artype = float
    M = max([N,D])
    Np = M
    Dp = M
    bprime = num.zeros(Np+1,artype)
    aprime = num.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] = num.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] = num.real(val)
    #return aprime, bprime
    return normalize(bprime, aprime)

 def centroid(self):
     '''Find the centroid x and y from these coefficients'''
     
     xcen = 0.0
     ycen = 0.0
     
     fluxtot = self.total_flux()
     
     for i1 in range(self.n1):            
         if i1%2==0: continue #consider odd i1
                         
         for i2 in range(self.n2):                    
             if i2%2!=0: continue # consider even i2
                                                                     
             xcen = xcen + np.power(i1+1,0.5)*np.power(2,0.5*(2-i1-i2))* np.power(comb(i1+1,(i1+1)/2)*comb(i2,i2/2),0.5)*self.coeff[i1,i2]
                 
             
     for i1 in range(self.n1):
         if i1%2!=0: continue #consider even i1   
                  
         for i2 in range(self.n2):
             if i2%2==0:continue # consider odd i2
                                                      
             ycen = ycen + np.power(i2+1,0.5)*np.power(2,0.5*(2-i2-i1))* np.power(comb(i2+1,(i2+1)/2)*comb(i1,i1/2),0.5)*self.coeff[i1,i2]
         
     xcen = xcen*np.sqrt(pi)*self.beta*self.beta/fluxtot                    
     ycen = ycen*np.sqrt(pi)*self.beta*self.beta/fluxtot
             
     return xcen,ycen

def multiple_alignment(word_list):
    '''Returns the multiple alignment of a given list of words.'''
    from itertools import product
    from operator import add, mul
    from scipy.misc import comb

    # There are some issues scoring the first symbols, so force a match here and remove it from the alignment later.
    word_list = ['$'+word for word in word_list]

    # Initialize scoring and backtrack dictionaries, along with the indices and base score.
    S, backtrack = {}, {}
    perm_list = list(product([0, -1], repeat=len(word_list)))[1:]
    base_score = -1*comb(len(word_list), 2, exact=True)

    for index in product(*map(xrange,map(lambda s: len(s) + 1, word_list))):

        # We forced a match with the first symbols, so the zero-shell should lead to the zero index.
        if reduce(mul, index) == 0:
            # Since we forced a match with the first symbol, we want to force starting point to be the zero index.
            if sum(index) == 0:
                # All symbols match.
                S[index] = 0
            else:
                # Make it smaller than the lowest possible score.
                S[index] = 2*base_score*reduce(add, map(len, word_list))

        else:
            # Use previous scores to determine the best score for the current index.
            previous_scores = [S[tuple(map(add, index, perm))] for perm in perm_list]
            current_index_scores = []
            for perm in perm_list:
                chars = [word_list[i][index[i]-1] if perm_value == -1 else '-' for i, perm_value in enumerate(perm)]
                current_index_scores.append(base_score + sum([comb(chars.count(ch), 2, exact=True) for ch in set(chars)]))

            scores = map(add, previous_scores, current_index_scores)
            backtrack[index], S[index] = max(enumerate(scores), key=lambda p: p[1])

    # Initialize the alignment and indicies.
    alignment = word_list
    current_index = map(len, word_list)

    # Get the max score.
    # Note: The forced match at start of each word does not change the max score, as matched symbols have a score of zero.
    max_score = S[tuple(current_index)]

    # Quick lambda function to insert indels.
    insert_indel = lambda word, i: word[:i] + '-' + word[i:]

    # Insert indels to get the alignment.
    while reduce(mul, current_index) != 0:
        for i, perm_value in enumerate(perm_list[backtrack[tuple(current_index)]]):
            if perm_value == 0:
                alignment[i] = insert_indel(alignment[i], current_index[i])
            else:
                current_index[i] -= 1

    # Note: We don't need to prepend any indels because we forced a match at the start of all words.
    # Remove the forced match from all alignments to recover the correct alignment.
    return [str(max_score)] + [str(aligned[1:].seq) for aligned in alignment]

def elm_comb(rho, N, m):
    offered_traffic = rho / (1 + rho)
    numerator = comb(N-1, m, exact=True)*pow(offered_traffic, m)*pow(1-offered_traffic, N-m)
    denominator = 0
    for i in range(0, m+1):
        denominator += comb(N-1, i, exact=True)*pow(offered_traffic,  i)*pow(1-offered_traffic, N-i)
    
    return numerator / denominator

 def rectangles_numbers(M, N):
     '''
     正方形是屬於長方形的一種。或者說，正方形是長方形的特例。矩形的定義是四個角皆為 90° 的四邊形。
     如果要算不包含正方形的長方形，只需把所有的算出，然後減去正方形的個數。
     '''
     horizontal = sm.comb(M+1, 2, exact=True)
     vertical = sm.comb(N+1, 2, exact=True)
     return horizontal * vertical

def hyperg(k, N, m, n, verbose=False):
    """
    k = intersection; m = population 1; n = population 2; N = total population.
    """
    exact=1
    if verbose and k%100==0: printnow("k: %s" % k)
    #return comb(m,k,exact) * comb(N-m,n-k,exact)/ comb(N,n,exact)
    return comb(m,k) * comb(N-m,n-k)/ comb(N,n)

 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 hypergeometricown(self, N, K, n, k):
   """
     N= total number of genes in population
     K= number of GOA
     n= select a sample (top 50, bottom half, etc.)
     k= number of successes in the sample
   """
   return comb(K, k) * comb(N-K, n-k) / comb(N, n)

def mv_hypergeometric(x, m):
    """
    x : number of draws for each category.
    m : size of each category.
    """
    x = np.asarray(x)
    m = np.asarray(m)
    return log(comb(m, x).prod() / comb(m.sum(), x.sum()))