treetime/treeregression.py
def find_best_root(self, force_positive=True, slope=None):
        """
        determine the position on the tree that minimizes the bilinear
        product of the inverse covariance and the data vectors.

        Returns
        -------
         best_root : (dict)
            dictionary with the node, the fraction `x` at which the branch
            is to be split, and the regression parameters
        """
        self._calculate_averages()
        best_root = {"chisq": np.inf}
        for n in self.tree.find_clades():
            if n==self.tree.root:
                continue

            tv = self.tip_value(n)
            bv = self.branch_value(n)
            var = self.branch_variance(n)
            x, chisq = self._optimal_root_along_branch(n, tv, bv, var, slope=slope)

            if (chisq<best_root["chisq"]):
                tmpQ = self.propagate_averages(n, tv, bv*x, var*x) \
                     + self.propagate_averages(n, tv, bv*(1-x), var*(1-x), outgroup=True)
                reg = base_regression(tmpQ, slope=slope)
                if reg["slope"]>=0 or (force_positive==False):
                    best_root = {"node":n, "split":x}
                    best_root.update(reg)

        if 'node' not in best_root:
            print("TreeRegression.find_best_root: No valid root found!", force_positive)
            return None

        if 'hessian' in best_root:
            # calculate differentials with respect to x
            deriv = []
            n = best_root["node"]
            tv = self.tip_value(n)
            bv = self.branch_value(n)
            var = self.branch_variance(n)
            for dx in [-0.001, 0.001]:
                y = min(1.0, max(0.0, best_root["split"]+dx))
                tmpQ = self.propagate_averages(n, tv, bv*y, var*y) \
                     + self.propagate_averages(n, tv, bv*(1-y), var*(1-y), outgroup=True)
                reg = base_regression(tmpQ, slope=slope)
                deriv.append([y,reg['chisq'], tmpQ[tavgii], tmpQ[davgii]])

            estimator_hessian = np.zeros((3,3))
            estimator_hessian[:2,:2] = best_root['hessian']
            estimator_hessian[2,2] = (deriv[0][1] + deriv[1][1] - 2.0*best_root['chisq'])/(deriv[0][0] - deriv[1][0])**2
            # estimator_hessian[2,0] = (deriv[0][2] - deriv[1][2])/(deriv[0][0] - deriv[1][0])
            # estimator_hessian[2,1] = (deriv[0][3] - deriv[1][3])/(deriv[0][0] - deriv[1][0])
            estimator_hessian[0,2] = estimator_hessian[2,0]
            estimator_hessian[1,2] = estimator_hessian[2,1]
            best_root['hessian'] = estimator_hessian
            best_root['cov'] = np.linalg.inv(estimator_hessian)

        return best_root
def find_best_root(self, force_positive=True, slope=None):
        self._calculate_averages()
        best_root = {"chisq": np.inf}
        for n in self.tree.find_clades():
            if n==self.tree.root:
                continue

            tv = self.tip_value(n)
            bv = self.branch_value(n)
            var = self.branch_variance(n)
            x, chisq = self._optimal_root_along_branch(n, tv, bv, var, slope=slope)

            if (chisq<best_root["chisq"]):
                tmpQ = self.propagate_averages(n, tv, bv*x, var*x) \
                     + self.propagate_averages(n, tv, bv*(1-x), var*(1-x), outgroup=True)
                reg = base_regression(tmpQ, slope=slope)
                if reg["slope"]>=0 or (force_positive==False):
                    best_root = {"node":n, "split":x}
                    best_root.update(reg)

        if 'node' not in best_root:
            print("TreeRegression.find_best_root: No valid root found!", force_positive)
            return None

        if 'hessian' in best_root:
            # calculate differentials with respect to x
            deriv = []
            n = best_root["node"]
            tv = self.tip_value(n)
            bv = self.branch_value(n)
            var = self.branch_variance(n)
            for dx in [-0.001, 0.001]:
                y = min(1.0, max(0.0, best_root["split"]+dx))
                tmpQ = self.propagate_averages(n, tv, bv*y, var*y) \
                     + self.propagate_averages(n, tv, bv*(1-y), var*(1-y), outgroup=True)
                reg = base_regression(tmpQ, slope=slope)
                deriv.append([y,reg['chisq'], tmpQ[tavgii], tmpQ[davgii]])

            estimator_hessian = np.zeros((3,3))
            estimator_hessian[:2,:2] = best_root['hessian']
            estimator_hessian[2,2] = (deriv[0][1] + deriv[1][1] - 2.0*best_root['chisq'])/(deriv[0][0] - deriv[1][0])**2
            # estimator_hessian[2,0] = (deriv[0][2] - deriv[1][2])/(deriv[0][0] - deriv[1][0])
            # estimator_hessian[2,1] = (deriv[0][3] - deriv[1][3])/(deriv[0][0] - deriv[1][0])
            estimator_hessian[0,2] = estimator_hessian[2,0]
            estimator_hessian[1,2] = estimator_hessian[2,1]
            best_root['hessian'] = estimator_hessian
            best_root['cov'] = np.linalg.inv(estimator_hessian)

        return best_root
def find_best_root ( self , force_positive = True , slope = None ) : self . _calculate_averages ( ) best_root = { "chisq" : np . inf } for n in self . tree . find_clades ( ) : if n == self . tree . root : continue tv = self . tip_value ( n ) bv = self . branch_value ( n ) var = self . branch_variance ( n ) x , chisq = self . _optimal_root_along_branch ( n , tv , bv , var , slope = slope ) if ( chisq < best_root [ "chisq" ] ) : tmpQ = self . propagate_averages ( n , tv , bv * x , var * x ) + self . propagate_averages ( n , tv , bv * ( 1 - x ) , var * ( 1 - x ) , outgroup = True ) reg = base_regression ( tmpQ , slope = slope ) if reg [ "slope" ] >= 0 or ( force_positive == False ) : best_root = { "node" : n , "split" : x } best_root . update ( reg ) if 'node' not in best_root : print ( "TreeRegression.find_best_root: No valid root found!" , force_positive ) return None if 'hessian' in best_root : # calculate differentials with respect to x deriv = [ ] n = best_root [ "node" ] tv = self . tip_value ( n ) bv = self . branch_value ( n ) var = self . branch_variance ( n ) for dx in [ - 0.001 , 0.001 ] : y = min ( 1.0 , max ( 0.0 , best_root [ "split" ] + dx ) ) tmpQ = self . propagate_averages ( n , tv , bv * y , var * y ) + self . propagate_averages ( n , tv , bv * ( 1 - y ) , var * ( 1 - y ) , outgroup = True ) reg = base_regression ( tmpQ , slope = slope ) deriv . append ( [ y , reg [ 'chisq' ] , tmpQ [ tavgii ] , tmpQ [ davgii ] ] ) estimator_hessian = np . zeros ( ( 3 , 3 ) ) estimator_hessian [ : 2 , : 2 ] = best_root [ 'hessian' ] estimator_hessian [ 2 , 2 ] = ( deriv [ 0 ] [ 1 ] + deriv [ 1 ] [ 1 ] - 2.0 * best_root [ 'chisq' ] ) / ( deriv [ 0 ] [ 0 ] - deriv [ 1 ] [ 0 ] ) ** 2 # estimator_hessian[2,0] = (deriv[0][2] - deriv[1][2])/(deriv[0][0] - deriv[1][0]) # estimator_hessian[2,1] = (deriv[0][3] - deriv[1][3])/(deriv[0][0] - deriv[1][0]) estimator_hessian [ 0 , 2 ] = estimator_hessian [ 2 , 0 ] estimator_hessian [ 1 , 2 ] = estimator_hessian [ 2 , 1 ] best_root [ 'hessian' ] = estimator_hessian best_root [ 'cov' ] = np . linalg . inv ( estimator_hessian ) return best_root
determine the position on the tree that minimizes the bilinear
        product of the inverse covariance and the data vectors.

        Returns
        -------
         best_root : (dict)
            dictionary with the node, the fraction `x` at which the branch
            is to be split, and the regression parameters
determine the position on the tree that minimizes the bilinear product of the inverse covariance and the data vectors .