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):
        """
        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
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