treetime/gtr.py
def optimal_t_compressed(self, seq_pair, multiplicity, profiles=False, tol=1e-10):
        """
        Find the optimal distance between the two sequences, for compressed sequences

        Parameters
        ----------

         seq_pair : compressed_sequence_pair
            Compressed representation of sequences along a branch, either
            as tuple of state pairs or as tuple of profiles.

         multiplicity : array
            Number of times each state pair in seq_pair appears (if profile==False)

            Number of times an alignment pattern is observed (if profiles==True)

         profiles : bool, default False
            The standard branch length optimization assumes fixed sequences at
            either end of the branch. With profiles==True, optimization is performed
            while summing over all possible states of the nodes at either end of the
            branch. Note that the meaning/format of seq_pair and multiplicity
            depend on the value of profiles.

        """

        def _neg_prob(t, seq_pair, multiplicity):
            """
            Probability to observe a child given the the parent state, transition
            matrix, and the time of evolution (branch length).

            Parameters
            ----------

             t : double
                Branch length (time between sequences)

             parent :  numpy.array
                Parent sequence

             child : numpy.array
                Child sequence

             tm :  GTR
                Model of evolution

            Returns
            -------

             prob : double
                Negative probability of the two given sequences
                to be separated by the time t.
            """
            if profiles:
                res = -1.0*self.prob_t_profiles(seq_pair, multiplicity,t**2, return_log=True)
                return res
            else:
                return -1.0*self.prob_t_compressed(seq_pair, multiplicity,t**2, return_log=True)

        try:
            from scipy.optimize import minimize_scalar
            opt = minimize_scalar(_neg_prob,
                    bounds=[-np.sqrt(ttconf.MAX_BRANCH_LENGTH),np.sqrt(ttconf.MAX_BRANCH_LENGTH)],
                    args=(seq_pair, multiplicity), tol=tol)
            new_len = opt["x"]**2
            if 'success' not in opt:
                opt['success'] = True
                self.logger("WARNING: the optimization result does not contain a 'success' flag:"+str(opt),4, warn=True)
        except:
            import scipy
            print('legacy scipy', scipy.__version__)
            from scipy.optimize import fminbound
            new_len = fminbound(_neg_prob,
                      -np.sqrt(ttconf.MAX_BRANCH_LENGTH),np.sqrt(ttconf.MAX_BRANCH_LENGTH),
                       args=(seq_pair, multiplicity))
            new_len = new_len**2
            opt={'success':True}

        if new_len > .9 * ttconf.MAX_BRANCH_LENGTH:
            self.logger("WARNING: GTR.optimal_t_compressed -- The branch length seems to be very long!", 4, warn=True)

        if opt["success"] != True:
            # return hamming distance: number of state pairs where state differs/all pairs
            new_len =  np.sum(multiplicity[seq_pair[:,1]!=seq_pair[:,0]])/np.sum(multiplicity)

        return new_len
def optimal_t_compressed(self, seq_pair, multiplicity, profiles=False, tol=1e-10):
        def _neg_prob(t, seq_pair, multiplicity):
            """
            Probability to observe a child given the the parent state, transition
            matrix, and the time of evolution (branch length).

            Parameters
            ----------

             t : double
                Branch length (time between sequences)

             parent :  numpy.array
                Parent sequence

             child : numpy.array
                Child sequence

             tm :  GTR
                Model of evolution

            Returns
            -------

             prob : double
                Negative probability of the two given sequences
                to be separated by the time t.
            """
            if profiles:
                res = -1.0*self.prob_t_profiles(seq_pair, multiplicity,t**2, return_log=True)
                return res
            else:
                return -1.0*self.prob_t_compressed(seq_pair, multiplicity,t**2, return_log=True)

        try:
            from scipy.optimize import minimize_scalar
            opt = minimize_scalar(_neg_prob,
                    bounds=[-np.sqrt(ttconf.MAX_BRANCH_LENGTH),np.sqrt(ttconf.MAX_BRANCH_LENGTH)],
                    args=(seq_pair, multiplicity), tol=tol)
            new_len = opt["x"]**2
            if 'success' not in opt:
                opt['success'] = True
                self.logger("WARNING: the optimization result does not contain a 'success' flag:"+str(opt),4, warn=True)
        except:
            import scipy
            print('legacy scipy', scipy.__version__)
            from scipy.optimize import fminbound
            new_len = fminbound(_neg_prob,
                      -np.sqrt(ttconf.MAX_BRANCH_LENGTH),np.sqrt(ttconf.MAX_BRANCH_LENGTH),
                       args=(seq_pair, multiplicity))
            new_len = new_len**2
            opt={'success':True}

        if new_len > .9 * ttconf.MAX_BRANCH_LENGTH:
            self.logger("WARNING: GTR.optimal_t_compressed -- The branch length seems to be very long!", 4, warn=True)

        if opt["success"] != True:
            # return hamming distance: number of state pairs where state differs/all pairs
            new_len =  np.sum(multiplicity[seq_pair[:,1]!=seq_pair[:,0]])/np.sum(multiplicity)

        return new_len
def optimal_t_compressed ( self , seq_pair , multiplicity , profiles = False , tol = 1e-10 ) : def _neg_prob ( t , seq_pair , multiplicity ) : """
            Probability to observe a child given the the parent state, transition
            matrix, and the time of evolution (branch length).

            Parameters
            ----------

             t : double
                Branch length (time between sequences)

             parent :  numpy.array
                Parent sequence

             child : numpy.array
                Child sequence

             tm :  GTR
                Model of evolution

            Returns
            -------

             prob : double
                Negative probability of the two given sequences
                to be separated by the time t.
            """ if profiles : res = - 1.0 * self . prob_t_profiles ( seq_pair , multiplicity , t ** 2 , return_log = True ) return res else : return - 1.0 * self . prob_t_compressed ( seq_pair , multiplicity , t ** 2 , return_log = True ) try : from scipy . optimize import minimize_scalar opt = minimize_scalar ( _neg_prob , bounds = [ - np . sqrt ( ttconf . MAX_BRANCH_LENGTH ) , np . sqrt ( ttconf . MAX_BRANCH_LENGTH ) ] , args = ( seq_pair , multiplicity ) , tol = tol ) new_len = opt [ "x" ] ** 2 if 'success' not in opt : opt [ 'success' ] = True self . logger ( "WARNING: the optimization result does not contain a 'success' flag:" + str ( opt ) , 4 , warn = True ) except : import scipy print ( 'legacy scipy' , scipy . __version__ ) from scipy . optimize import fminbound new_len = fminbound ( _neg_prob , - np . sqrt ( ttconf . MAX_BRANCH_LENGTH ) , np . sqrt ( ttconf . MAX_BRANCH_LENGTH ) , args = ( seq_pair , multiplicity ) ) new_len = new_len ** 2 opt = { 'success' : True } if new_len > .9 * ttconf . MAX_BRANCH_LENGTH : self . logger ( "WARNING: GTR.optimal_t_compressed -- The branch length seems to be very long!" , 4 , warn = True ) if opt [ "success" ] != True : # return hamming distance: number of state pairs where state differs/all pairs new_len = np . sum ( multiplicity [ seq_pair [ : , 1 ] != seq_pair [ : , 0 ] ] ) / np . sum ( multiplicity ) return new_len
Find the optimal distance between the two sequences, for compressed sequences

        Parameters
        ----------

         seq_pair : compressed_sequence_pair
            Compressed representation of sequences along a branch, either
            as tuple of state pairs or as tuple of profiles.

         multiplicity : array
            Number of times each state pair in seq_pair appears (if profile==False)

            Number of times an alignment pattern is observed (if profiles==True)

         profiles : bool, default False
            The standard branch length optimization assumes fixed sequences at
            either end of the branch. With profiles==True, optimization is performed
            while summing over all possible states of the nodes at either end of the
            branch. Note that the meaning/format of seq_pair and multiplicity
            depend on the value of profiles.
Find the optimal distance between the two sequences for compressed sequences