treetime/gtr_site_specific.py
def infer(cls, sub_ija, T_ia, root_state, pc=0.01,
              gap_limit=0.01, Nit=30, dp=1e-5, **kwargs):
        """
        Infer a GTR model by specifying the number of transitions and time spent in each
        character. The basic equation that is being solved is

        :math:`n_{ij} = pi_i W_{ij} T_j`

        where :math:`n_{ij}` are the transitions, :math:`pi_i` are the equilibrium
        state frequencies, :math:`W_{ij}` is the "substitution attempt matrix",
        while :math:`T_i` is the time on the tree spent in character state
        :math:`i`. To regularize the process, we add pseudocounts and also need
        to account for the fact that the root of the tree is in a particular
        state. the modified equation is

        :math:`n_{ij} + pc = pi_i W_{ij} (T_j+pc+root\_state)`

        Parameters
        ----------

         nija : nxn matrix
            The number of times a change in character state is observed
            between state j and i at position a

         Tia :n vector
            The time spent in each character state at position a

         root_state : n vector
            The number of characters in state i in the sequence
            of the root node.

         pc : float
            Pseudocounts, this determines the lower cutoff on the rate when
            no substitutions are observed

         **kwargs:
            Key word arguments to be passed

        Keyword Args
        ------------

         alphabet : str
            Specify alphabet when applicable. If the alphabet specification
            is required, but no alphabet is specified, the nucleotide alphabet will be used as default.

        """
        from scipy import linalg as LA
        gtr = cls(**kwargs)
        gtr.logger("GTR: model inference ",1)
        q = len(gtr.alphabet)
        L = sub_ija.shape[-1]

        n_iter = 0
        n_ija = np.copy(sub_ija)
        n_ija[range(q),range(q),:] = 0
        n_ij = n_ija.sum(axis=-1)

        m_ia = np.sum(n_ija,axis=1) + root_state + pc
        n_a = n_ija.sum(axis=1).sum(axis=0) + pc

        Lambda = np.sum(root_state,axis=0) + q*pc
        p_ia_old=np.zeros((q,L))
        p_ia = np.ones((q,L))/q
        mu_a = np.ones(L)

        W_ij = np.ones((q,q)) - np.eye(q)

        while (LA.norm(p_ia_old-p_ia)>dp) and n_iter<Nit:
            n_iter += 1
            p_ia_old = np.copy(p_ia)
            S_ij = np.einsum('a,ia,ja',mu_a, p_ia, T_ia)
            W_ij = (n_ij + n_ij.T + pc)/(S_ij + S_ij.T + pc)

            avg_pi = p_ia.mean(axis=-1)
            average_rate = W_ij.dot(avg_pi).dot(avg_pi)
            W_ij = W_ij/average_rate
            mu_a *=average_rate

            p_ia = m_ia/(mu_a*np.dot(W_ij,T_ia)+Lambda)
            p_ia = p_ia/p_ia.sum(axis=0)

            mu_a = n_a/(pc+np.einsum('ia,ij,ja->a', p_ia, W_ij, T_ia))

        if n_iter >= Nit:
            gtr.logger('WARNING: maximum number of iterations has been reached in GTR inference',3, warn=True)
            if LA.norm(p_ia_old-p_ia) > dp:
                gtr.logger('the iterative scheme has not converged',3,warn=True)
        if gtr.gap_index is not None:
            for p in range(p_ia.shape[-1]):
                if p_ia[gtr.gap_index,p]<gap_limit:
                    gtr.logger('The model allows for gaps which are estimated to occur at a low fraction of %1.3e'%p_ia[gtr.gap_index]+
                           '\n\t\tthis can potentially result in artificats.'+
                           '\n\t\tgap fraction will be set to %1.4f'%gap_limit,2,warn=True)
                p_ia[gtr.gap_index,p] = gap_limit
                p_ia[:,p] /= p_ia[:,p].sum()

        gtr.assign_rates(mu=mu_a, W=W_ij, pi=p_ia)
        return gtr
def infer(cls, sub_ija, T_ia, root_state, pc=0.01,
              gap_limit=0.01, Nit=30, dp=1e-5, **kwargs):
        """
        Infer a GTR model by specifying the number of transitions and time spent in each
        character. The basic equation that is being solved is

        :math:`n_{ij} = pi_i W_{ij} T_j`

        where :math:`n_{ij}` are the transitions, :math:`pi_i` are the equilibrium
        state frequencies, :math:`W_{ij}` is the "substitution attempt matrix",
        while :math:`T_i` is the time on the tree spent in character state
        :math:`i`. To regularize the process, we add pseudocounts and also need
        to account for the fact that the root of the tree is in a particular
        state. the modified equation is

        :math:`n_{ij} + pc = pi_i W_{ij} (T_j+pc+root\_state)`

        Parameters
        ----------

         nija : nxn matrix
            The number of times a change in character state is observed
            between state j and i at position a

         Tia :n vector
            The time spent in each character state at position a

         root_state : n vector
            The number of characters in state i in the sequence
            of the root node.

         pc : float
            Pseudocounts, this determines the lower cutoff on the rate when
            no substitutions are observed

         **kwargs:
            Key word arguments to be passed

        Keyword Args
        ------------

         alphabet : str
            Specify alphabet when applicable. If the alphabet specification
            is required, but no alphabet is specified, the nucleotide alphabet will be used as default.

        """
        from scipy import linalg as LA
        gtr = cls(**kwargs)
        gtr.logger("GTR: model inference ",1)
        q = len(gtr.alphabet)
        L = sub_ija.shape[-1]

        n_iter = 0
        n_ija = np.copy(sub_ija)
        n_ija[range(q),range(q),:] = 0
        n_ij = n_ija.sum(axis=-1)

        m_ia = np.sum(n_ija,axis=1) + root_state + pc
        n_a = n_ija.sum(axis=1).sum(axis=0) + pc

        Lambda = np.sum(root_state,axis=0) + q*pc
        p_ia_old=np.zeros((q,L))
        p_ia = np.ones((q,L))/q
        mu_a = np.ones(L)

        W_ij = np.ones((q,q)) - np.eye(q)

        while (LA.norm(p_ia_old-p_ia)>dp) and n_iter<Nit:
            n_iter += 1
            p_ia_old = np.copy(p_ia)
            S_ij = np.einsum('a,ia,ja',mu_a, p_ia, T_ia)
            W_ij = (n_ij + n_ij.T + pc)/(S_ij + S_ij.T + pc)

            avg_pi = p_ia.mean(axis=-1)
            average_rate = W_ij.dot(avg_pi).dot(avg_pi)
            W_ij = W_ij/average_rate
            mu_a *=average_rate

            p_ia = m_ia/(mu_a*np.dot(W_ij,T_ia)+Lambda)
            p_ia = p_ia/p_ia.sum(axis=0)

            mu_a = n_a/(pc+np.einsum('ia,ij,ja->a', p_ia, W_ij, T_ia))

        if n_iter >= Nit:
            gtr.logger('WARNING: maximum number of iterations has been reached in GTR inference',3, warn=True)
            if LA.norm(p_ia_old-p_ia) > dp:
                gtr.logger('the iterative scheme has not converged',3,warn=True)
        if gtr.gap_index is not None:
            for p in range(p_ia.shape[-1]):
                if p_ia[gtr.gap_index,p]<gap_limit:
                    gtr.logger('The model allows for gaps which are estimated to occur at a low fraction of %1.3e'%p_ia[gtr.gap_index]+
                           '\n\t\tthis can potentially result in artificats.'+
                           '\n\t\tgap fraction will be set to %1.4f'%gap_limit,2,warn=True)
                p_ia[gtr.gap_index,p] = gap_limit
                p_ia[:,p] /= p_ia[:,p].sum()

        gtr.assign_rates(mu=mu_a, W=W_ij, pi=p_ia)
        return gtr
def infer ( cls , sub_ija , T_ia , root_state , pc = 0.01 , gap_limit = 0.01 , Nit = 30 , dp = 1e-5 , * * kwargs ) : from scipy import linalg as LA gtr = cls ( * * kwargs ) gtr . logger ( "GTR: model inference " , 1 ) q = len ( gtr . alphabet ) L = sub_ija . shape [ - 1 ] n_iter = 0 n_ija = np . copy ( sub_ija ) n_ija [ range ( q ) , range ( q ) , : ] = 0 n_ij = n_ija . sum ( axis = - 1 ) m_ia = np . sum ( n_ija , axis = 1 ) + root_state + pc n_a = n_ija . sum ( axis = 1 ) . sum ( axis = 0 ) + pc Lambda = np . sum ( root_state , axis = 0 ) + q * pc p_ia_old = np . zeros ( ( q , L ) ) p_ia = np . ones ( ( q , L ) ) / q mu_a = np . ones ( L ) W_ij = np . ones ( ( q , q ) ) - np . eye ( q ) while ( LA . norm ( p_ia_old - p_ia ) > dp ) and n_iter < Nit : n_iter += 1 p_ia_old = np . copy ( p_ia ) S_ij = np . einsum ( 'a,ia,ja' , mu_a , p_ia , T_ia ) W_ij = ( n_ij + n_ij . T + pc ) / ( S_ij + S_ij . T + pc ) avg_pi = p_ia . mean ( axis = - 1 ) average_rate = W_ij . dot ( avg_pi ) . dot ( avg_pi ) W_ij = W_ij / average_rate mu_a *= average_rate p_ia = m_ia / ( mu_a * np . dot ( W_ij , T_ia ) + Lambda ) p_ia = p_ia / p_ia . sum ( axis = 0 ) mu_a = n_a / ( pc + np . einsum ( 'ia,ij,ja->a' , p_ia , W_ij , T_ia ) ) if n_iter >= Nit : gtr . logger ( 'WARNING: maximum number of iterations has been reached in GTR inference' , 3 , warn = True ) if LA . norm ( p_ia_old - p_ia ) > dp : gtr . logger ( 'the iterative scheme has not converged' , 3 , warn = True ) if gtr . gap_index is not None : for p in range ( p_ia . shape [ - 1 ] ) : if p_ia [ gtr . gap_index , p ] < gap_limit : gtr . logger ( 'The model allows for gaps which are estimated to occur at a low fraction of %1.3e' % p_ia [ gtr . gap_index ] + '\n\t\tthis can potentially result in artificats.' + '\n\t\tgap fraction will be set to %1.4f' % gap_limit , 2 , warn = True ) p_ia [ gtr . gap_index , p ] = gap_limit p_ia [ : , p ] /= p_ia [ : , p ] . sum ( ) gtr . assign_rates ( mu = mu_a , W = W_ij , pi = p_ia ) return gtr
Infer a GTR model by specifying the number of transitions and time spent in each
        character. The basic equation that is being solved is

        :math:`n_{ij} = pi_i W_{ij} T_j`

        where :math:`n_{ij}` are the transitions, :math:`pi_i` are the equilibrium
        state frequencies, :math:`W_{ij}` is the "substitution attempt matrix",
        while :math:`T_i` is the time on the tree spent in character state
        :math:`i`. To regularize the process, we add pseudocounts and also need
        to account for the fact that the root of the tree is in a particular
        state. the modified equation is

        :math:`n_{ij} + pc = pi_i W_{ij} (T_j+pc+root\_state)`

        Parameters
        ----------

         nija : nxn matrix
            The number of times a change in character state is observed
            between state j and i at position a

         Tia :n vector
            The time spent in each character state at position a

         root_state : n vector
            The number of characters in state i in the sequence
            of the root node.

         pc : float
            Pseudocounts, this determines the lower cutoff on the rate when
            no substitutions are observed

         **kwargs:
            Key word arguments to be passed

        Keyword Args
        ------------

         alphabet : str
            Specify alphabet when applicable. If the alphabet specification
            is required, but no alphabet is specified, the nucleotide alphabet will be used as default.
Infer a GTR model by specifying the number of transitions and time spent in each character . The basic equation that is being solved is