def grow(self):
        """Grow the population to carrying capacity
        
        The final population size is determined based on the proportion of
        producers present. This population is determined by drawing from a
        multinomial with the probability of each genotype proportional to its
        abundance times its fitness.
        """

        if self.is_empty():
            return

        if not self.diluted:
            return

        landscape = self.metapopulation.fitness_landscape

        final_size = self.capacity_min + \
                (self.capacity_max - self.capacity_min) * \
                self.prop_producers()

        grow_probs = self.abundances * (landscape/nsum(landscape))

        if nsum(grow_probs) > 0:
            norm_grow_probs = grow_probs/nsum(grow_probs)
            self.abundances = multinomial(final_size, norm_grow_probs, 1)[0]

def get_pos_norm(pos):
    """ Translate positions so that centroid is in the origin and return mean
    norm of the translated positions. """
    pos = Positions.create(pos)
    assert(len(pos) == 1)
    n = len(pos[0])
    p = zeros((n, 2))
    for i, node in enumerate(pos[0]):
        p[i, :] = pos[0][node]
    centroid = p.sum(axis=0)/n
    p -= tile(centroid, (n, 1))
    p_norm = nsum(sqrt(nsum(square(p), axis=1)))/n
    return p_norm

 def __init__(self, hamiltonian, beta, mu, verbose = False):
     self.mu = mu
     CanonicalEnsemble.__init__(self, hamiltonian, beta, verbose)
     c = AnnihilationOperator(self.singleParticleBasis)
     muMatrix = mu * nsum([c[orb].H.dot(c[orb]) for orb in self.orderedSingleParticleStates], axis = 0)
     self.hamiltonian.matrix = self.hamiltonian.matrix - muMatrix
     self.filling = None

    def RHS(
        self, N, t
    ):
        dNdt = zeros_like(N)

        if self.gamma is not None and self.betadxi is not None:
            # Death breakup term
            dNdt[1:] -= N[1:] * self.gamma[1:]
            dNdt[:-1] += self.nu * dot(
                self.betadxi[:-1, 1:], N[1:] * self.gamma[1:])

        if self.Q is not None:
            Cd = zeros_like(dNdt)
            for i in arange(self.number_of_classes / 2):
                ind = slice(i, self.number_of_classes - i - 1)
                Cb = self.Q[i, ind] * N[i] * N[ind]
                Cd[i] += nsum(Cb)
                Cd[(i + 1):(self.number_of_classes - i - 1)] += Cb[1:]
                Cb[0] = 0.5 * Cb[0]
                dNdt[(2 * i + 1):] += Cb

            dNdt -= Cd
        if self.theta is not None:
            dNdt += (self.n0 * self.A0 - N / self.theta)

        # print('Time = {0:g}'.format(t))
        return dNdt

def sig_q_e_LHS( e ):
    ''' Summation / integration  over the random domain '''
    q_e_grid = q( e, T_la[:, None, None, None, None], T_xi[None, :, None, None, None],
                   T_E[None, None, :, None, None], T_th[None, None, None, :, None],
                   T_A[None, None, None, None, :] )
    q_dG_grid = q_e_grid ** 2 / n_int ** n_k
    return sqrt( nsum( q_dG_grid ) - mu_q_e_LHS( e ) ** 2 )

def sig_q_e( e ):
    ''' Summation / integration  over the random domain '''
    q_e_grid = q( e, Theta_la[:, None, None, None, None], Theta_xi[None, :, None, None, None],
                   Theta_E[None, None, :, None, None], Theta_th[None, None, None, :, None],
                   Theta_A[None, None, None, None, :] )
    q_dG_grid = q_e_grid ** 2 * dG_grid
    return sqrt( nsum( q_dG_grid ) - mu_q_e( e ) ** 2 )

def mu_q_e( e ):
    ''' Summation / integration  over the random domain '''
    q_e_grid = q( e, Theta_la[:, None, None, None, None], Theta_xi[None, :, None, None, None],
                   Theta_E[None, None, :, None, None], Theta_th[None, None, None, :, None],
                   Theta_A[None, None, None, None, :] )
    q_dG_grid = q_e_grid * dG_grid # element by element product of two (m,m) arrays
    return nsum( q_dG_grid ) # nsum has been imported at line 3 from numpy 

 def __init__(self, blocksizes, all_real = True):
     Blocks.__init__(self, blocksizes)
     for size in blocksizes:
         if all_real:
             self.datablocks.append(asmatrix(zeros([size, size])))
         else:
             self.datablocks.append(asmatrix(zeros([size, size], dtype=complex)))
     self.shape = [nsum(blocksizes)]*2

 def likelihood_func_deriv(theta, N, X, B, z_A, z_B, mutation):
     # derivative of the likelihood of the function for theta
     p_A, p_B, f_A, f_B = p_read(theta, N, X, B, mutation)
     p_A_deriv = B * (1.0-B) * ((1-theta+B*theta) **(-2.0))
     p_B_deriv = ((1.0-B*theta) * (-B) - (B-B*theta) * (-B))/ ((1.0-B*theta)**2)
     f_A_deriv = nchoosek(N,X) * X *(p_A_deriv) * ((1.0-p_A)**(N-X)) + nchoosek(N,X) * (p_A**X) * (N-X) * (-p_A_deriv)
     f_B_deriv = nchoosek(N,X) * X *(p_B_deriv) * ((1.0-p_B)**(N-X)) + nchoosek(N,X) * (p_B**X) * (N-X) * (-p_B_deriv)
     l_deriv = -1.0*nsum(1.0/(z_A*f_A + z_B*f_B) * (z_A * f_A_deriv + z_B * f_B_deriv))
     return l_deriv

def clean_correlation_matrix(evals, evecs, max_eig, phylogeny=True):
    """ Cleans the correlation matrix of noise, and provides an option to remove the largest eigenvector to clean the matrix of phylogeny.

    Arguments:
    evals -- Eigenvalues of correlation matrix.
    evecs -- Eigenvectors of correlation matrix.
    max_eig -- Theoretical random maximum eigenvalue. We ignore anything below the minimum eigenvalue.
    """
    return real(nsum((x*outer(y,y) for x,y in zip(evals, evecs))))

def setHubbardMatrix(t, u, spins, orbitals, siteSpaceTransformation = None): # TODO rm siteSTrafo
    spins = range(len(spins))
    c = AnnihilationOperator([spins, orbitals])
    no = len(orbitals)
    ns = len(spins)
    spininds = range(len(spins))
    uMatrix = zeros([no,no,no,no,ns,ns])
    for i, j, k, l, s1, s2 in product(orbitals,orbitals,orbitals,orbitals,spins,spins):
        if i == k and j == l and i == j and s1 != s2:
            uMatrix[i, j, k, l, s1, s2] = u * .5
    if siteSpaceTransformation != None:
        p = array(siteSpaceTransformation)
        t = p.transpose().dot(t).dot(p)
        temp = uMatrix.copy()
        for i, j, k, l, s1, s2 in product(orbitals,orbitals,orbitals,orbitals,spins,spins):
            uMatrix[i,j,k,l,s1,s2] = nsum([p[i,m] * p[j,n] * temp[m,n,o,q,s1,s2] * p.transpose()[o,l] * p.transpose()[q,k] for m,n,o,q in product(orbitals,orbitals,orbitals,orbitals)], axis = 0)
    ht = [t[i,j] * c[s,i].H.dot(c[s,j]) for s in spins for i,j in product(orbitals, orbitals) if t[i,j] != 0]
    hu = [uMatrix[i, j, k, l, s1, s2] * c[s1,i].H.dot(c[s2, j].H).dot(c[s2, l]).dot(c[s1, k]) for i,j,k,l in product(orbitals, orbitals, orbitals, orbitals) for s1, s2 in product(spins, spins) if s1 != s2 and uMatrix[i, j, k, l, s1, s2] != 0]
    return nsum(ht + hu, axis = 0)

 def getFockspaceNr(self, occupationOfSingleParticleStates = None, singleParticleStateNr = None, singleParticleState = None):
     """Fock state number."""
     if occupationOfSingleParticleStates != None:
         return nsum([int(occ)*2**i for i, occ in enumerate(occupationOfSingleParticleStates)])
     elif singleParticleStateNr != None:
         return 2**singleParticleStateNr
     elif singleParticleState != None:
         return 2**self.getSingleParticleStateNr(*singleParticleState)
     else:
         assert False, 'Need parameter.'

 def setMu(self, mu):
     c = AnnihilationOperator(self.singleParticleBasis)
     nMatrix = nsum([c[orb].H.dot(c[orb]) for orb in self.orderedSingleParticleStates], axis = 0)
     self.hamiltonian.matrix = self.hamiltonian.matrix + self.mu * nMatrix
     self.mu = mu
     self.hamiltonian.matrix = self.hamiltonian.matrix - mu * nMatrix
     self.energyEigenvalues = None
     self.energyEigenstates = None
     self.partitionFunction = None
     self.occupation = dict()
     report('Chemical potential set to '+str(mu), self.verbose)

 def __getitem__(self, spState):
     """The single particle state is given by a tuple of quantum numbers. Returns scipy.sparse.coo_matrix"""
     instates = list()
     outstates = list()
     spStateOR = self.getOccupationRep(singleParticleState = spState)
     for fockStateNr in range(self.fockspaceSize):
         instateOR = self.getOccupationRep(fockStateNr)
         if instateOR[self.orderedSingleParticleStates.index(spState)] == '1':
             instates.append(fockStateNr)
             outstates.append(self.getFockspaceNr(annihilateOccRep(spStateOR, instateOR)))
     signs = [(-1)**nsum([1 for k in range(self.getSingleParticleStateNr(*spState)) if self.getOccupationRep(fockstateNr)[k] == '1']) for fockstateNr in instates]
     return coo_matrix((signs, (outstates, instates)), [self.fockspaceSize]*2)

    def __init__(
            self, number_of_classes, t, dxi, N0=None, xi0=None,
            beta=None, gamma=None, Q=None,
            theta=None, n0=None, A0=None):
        self.number_of_classes = number_of_classes
        if xi0 is None:
            self.xi0 = dxi
        else:
            self.xi0 = xi0
        self.n0 = n0
        self.theta = theta
        # Uniform grid
        self.xi = self.xi0 + dxi * arange(self.number_of_classes)
        if N0 is None:
            N0 = zeros_like(self.xi)
        else:
            N0 = array([
                quad(N0, self.xi[i] - dxi / 2., self.xi[i] + dxi / 2.)[0]
                for i in range(number_of_classes)])

        self.nu = 2.0  # Binary breakup
        # Kernels setup
        if gamma is not None:
            self.gamma = gamma(self.xi)
            self.betadxi = zeros(
                (self.number_of_classes, self.number_of_classes))
            for i in range(1, len(self.xi)):
                for j in range(i):
                    self.betadxi[j, i] = beta(self.xi[j], self.xi[i])
                self.betadxi[:, i] =\
                    self.betadxi[:, i] / nsum(self.betadxi[:, i])

        else:
            self.gamma = None
            self.betadxi = None

        if Q is not None:
            self.Q = zeros((self.number_of_classes, self.number_of_classes))
            for i in range(len(self.xi)):
                for j in range(len(self.xi)):
                    self.Q[i, j] = Q(self.xi[i], self.xi[j])  #
        else:
            self.Q = None

        if A0 is None:
            self.A0 = None
        else:
            self.A0 = array([
                quad(A0, self.xi[i] - dxi / 2., self.xi[i] + dxi / 2.)[0]
                for i in range(number_of_classes)])
        # Solve procedure
        self.N = odeint(lambda NN, t: self.RHS(NN, t), N0, t)

def lehmannSumDynamic(elementProducts, energyDifferences, partitionFunction, mesh, zeroFrequencyTerms, imaginaryOffset, nominatorCoefficients):
    if type(imaginaryOffset) != list:
        imaginaryOffset = [imaginaryOffset]*2
    result = dict()
    for statePair, nominators, denominators in zip(elementProducts.keys(), elementProducts.values(), energyDifferences.values()):
        data_p = list()
        for w in scatter_list(mesh):
            terms = [coeff * nom/(w + denom + complex(0, offset)) for noms, denom in zip(nominators, denominators) for nom, coeff, offset in zip(noms, nominatorCoefficients, imaginaryOffset)]
            if len(zeroFrequencyTerms.values()) > 0 and w == 0:
                terms += zeroFrequencyTerms[statePair]
            data_p.append(nsum(terms/partitionFunction))
        result.update({statePair: array(allgather_list(data_p))})
    return result

def do_pcoa(dists):
    'It does a Principal Coordinate Analysis on a distance matrix'
    # the code for this function is taken from pycogent metric_scaling.py
    # Principles of Multivariate analysis: A User's Perspective.
    # W.J. Krzanowski Oxford University Press, 2000. p106.

    dists = squareform(dists)

    e_matrix = (dists * dists) / -2.0
    f_matrix = _make_f_matrix(e_matrix)

    eigvals, eigvecs = eigh(f_matrix)
    eigvecs = eigvecs.transpose()
    # drop imaginary component, if we got one
    eigvals, eigvecs = eigvals.real, eigvecs.real

    # convert eigvals and eigvecs to point matrix
    # normalized eigenvectors with eigenvalues

    # get the coordinates of the n points on the jth axis of the Euclidean
    # representation as the elements of (sqrt(eigvalj))eigvecj
    # must take the absolute value of the eigvals since they can be negative
    pca_matrix = eigvecs * sqrt(abs(eigvals))[:, newaxis]

    # output
    # get order to output eigenvectors values. reports the eigvecs according
    # to their cooresponding eigvals from greatest to least
    vector_order = list(argsort(eigvals))
    vector_order.reverse()

    eigvals = eigvals[vector_order]

    # eigenvalues
    pcnts = (eigvals / nsum(eigvals)) * 100.0

    # the outputs
    # eigenvectors in the original pycogent implementation, here we name them
    # princoords
    # I think that we're doing: if the eigenvectors are written as columns,
    # the rows of the resulting table are the coordinates of the objects in
    # PCO space
    projections = []
    for name_i in range(dists.shape[0]):
        eigvect = [pca_matrix[vec_i, name_i] for vec_i in vector_order]
        projections.append(eigvect)
    projections = array(projections)

    return {'projections': projections,
            'var_percentages': pcnts}

 def setMuByFilling(self, filling, muMin, muMax, muTol = .001, maxiter = 100):
     c = AnnihilationOperator(self.singleParticleBasis)
     nMatrix = nsum([c[orb].H.dot(c[orb]) for orb in self.orderedSingleParticleStates], axis = 0)
     self.hamiltonian.matrix = self.hamiltonian.matrix + self.mu * nMatrix
     def fillingFunction(muTrial):
         self.hamiltonian.matrix = self.hamiltonian.matrix - muTrial * nMatrix
         self.calcEigensystem()
         self.calcPartitionFunction()
         self.calcOccupation()
         fillingTrial = self.getTotalOccupation()
         self.hamiltonian.matrix = self.hamiltonian.matrix + muTrial * nMatrix
         report('Filling(mu='+str(muTrial)+') = '+str(fillingTrial), self.verbose)
         self.filling = fillingTrial
         return fillingTrial - filling
     mu = bisect(fillingFunction, muMin, muMax, xtol = muTol, maxiter = maxiter)
     self.mu = mu
     self.hamiltonian.matrix = self.hamiltonian.matrix - mu * nMatrix
     report('Chemical potential set to '+str(mu), self.verbose)

def embedV(subspaceVectors, blocksizes):
    fockspaceSize = nsum(blocksizes)
    assert fockspaceSize == len(subspaceVectors), 'embedding will fail'
    iBlock = 0
    iBlockOrigin = 0
    vectors = zeros([len(subspaceVectors), fockspaceSize])
    x = list()
    y = list()
    data = list()
    for i, v in enumerate(subspaceVectors):
        for j, vj in enumerate(v):
            if not equals(vj, 0):
                y.append(j + iBlockOrigin) # TODO understand row/col exchange
                x.append(i)
                data.append(vj)
        if i == iBlockOrigin + blocksizes[iBlock] - 1:
            iBlockOrigin += blocksizes[iBlock]
            iBlock += 1
    return coo_matrix((data, (x,y)), [fockspaceSize]*2)

def _make_f_matrix(matrix):
    """It takes an E matrix and returns an F matrix

    The input is the output of make_E_matrix

    For each element in matrix subtract mean of corresponding row and
    column and add the mean of all elements in the matrix
    """
    num_rows, num_cols = matrix.shape
    # make a vector of the means for each row and column
    # column_means = (add.reduce(E_matrix) / num_rows)
    column_means = (add.reduce(matrix) / num_rows)[:, newaxis]
    trans_matrix = transpose(matrix)
    row_sums = add.reduce(trans_matrix)
    row_means = row_sums / num_cols
    # calculate the mean of the whole matrix
    matrix_mean = nsum(row_sums) / (num_rows * num_cols)
    # adjust each element in the E matrix to make the F matrix

    matrix -= row_means
    matrix -= column_means
    matrix += matrix_mean

    return matrix