treetime/treeregression.py
def clock_plot(self, add_internal=False, ax=None, regression=None,
                   confidence=True, n_sigma = 2, fs=14):
        """Plot root-to-tip distance vs time as a basic time-tree diagnostic

        Parameters
        ----------
        add_internal : bool, optional
            add internal nodes. this will only work if the tree has been dated already
        ax : None, optional
            an matplotlib axis to plot into. if non provided, a new figure is opened
        regression : None, optional
            a dict containing parameters of a root-to-tip vs time regression as
            returned by the function base_regression
        confidence : bool, optional
            add confidence area to the regression line
        n_sigma : int, optional
            number of standard deviations for the confidence area.
        fs : int, optional
            fontsize

        """
        import matplotlib.pyplot as plt
        if ax is None:
            plt.figure()
            ax=plt.subplot(111)

        self.tree.root._v=0
        for n in self.tree.get_nonterminals(order='preorder'):
            for c in n:
                c._v = n._v + self.branch_value(c)

        tips = self.tree.get_terminals()
        internal = self.tree.get_nonterminals()

        # get values of terminals
        xi = np.array([self.tip_value(n) for n in tips])
        yi = np.array([n._v for n in tips])
        ind = np.array([n.bad_branch  if hasattr(n, 'bad_branch') else False for n in tips])
        if add_internal:
            xi_int = np.array([n.numdate for n in internal])
            yi_int = np.array([n._v for n in internal])
            ind_int = np.array([n.bad_branch  if hasattr(n, 'bad_branch') else False  for n in internal])

        if regression:
            # plot regression line
            t_mrca = -regression['intercept']/regression['slope']
            if add_internal:
                time_span = np.max(xi_int[~ind_int]) - np.min(xi_int[~ind_int])
                x_vals = np.array([max(np.min(xi_int[~ind_int]), t_mrca) - 0.1*time_span, np.max(xi[~ind])+0.05*time_span])
            else:
                time_span = np.max(xi[~ind]) - np.min(xi[~ind])
                x_vals = np.array([max(np.min(xi[~ind]), t_mrca) - 0.1*time_span, np.max(xi[~ind]+0.05*time_span)])

            # plot confidence interval
            if confidence and 'cov' in regression:
                x_vals = np.linspace(x_vals[0], x_vals[1], 100)
                y_vals = regression['slope']*x_vals + regression['intercept']
                dev = n_sigma*np.array([np.sqrt(regression['cov'][:2,:2].dot(np.array([x, 1])).dot(np.array([x,1]))) for x in x_vals])
                dev_slope = n_sigma*np.sqrt(regression['cov'][0,0])
                ax.fill_between(x_vals, y_vals-dev, y_vals+dev, alpha=0.2)
                dp = np.array([regression['intercept']/regression['slope']**2,-1./regression['slope']])
                dev_rtt = n_sigma*np.sqrt(regression['cov'][:2,:2].dot(dp).dot(dp))

            else:
                dev_rtt = None
                dev_slope = None

            ax.plot(x_vals, regression['slope']*x_vals + regression['intercept'],
                    label = r"$y=\alpha + \beta t$"+"\n"+
                            r"$\beta=$%1.2e"%(regression["slope"])
                            + ("+/- %1.e"%dev_slope if dev_slope else "") +
                            "\nroot date: %1.1f"%(-regression['intercept']/regression['slope']) +
                            ("+/- %1.2f"%dev_rtt if dev_rtt else ""))


        ax.scatter(xi[~ind], yi[~ind], label=("tips" if add_internal else None))
        if ind.sum():
            try:
                # note: this is treetime specific
                tmp_x = np.array([np.mean(n.raw_date_constraint) if n.raw_date_constraint else None
                                  for n in self.tree.get_terminals()])
                ax.scatter(tmp_x[ind], yi[ind], label="ignored tips", c='r')
            except:
                pass
        if add_internal:
            ax.scatter(xi_int[~ind_int], yi_int[~ind_int], label="internal nodes")

        ax.set_ylabel('root-to-tip distance', fontsize=fs)
        ax.set_xlabel('date', fontsize=fs)
        ax.ticklabel_format(useOffset=False)
        ax.tick_params(labelsize=fs*0.8)
        ax.set_ylim([0, 1.1*np.max(yi)])
        plt.tight_layout()
        plt.legend(fontsize=fs*0.8)
def clock_plot(self, add_internal=False, ax=None, regression=None,
                   confidence=True, n_sigma = 2, fs=14):
        """Plot root-to-tip distance vs time as a basic time-tree diagnostic

        Parameters
        ----------
        add_internal : bool, optional
            add internal nodes. this will only work if the tree has been dated already
        ax : None, optional
            an matplotlib axis to plot into. if non provided, a new figure is opened
        regression : None, optional
            a dict containing parameters of a root-to-tip vs time regression as
            returned by the function base_regression
        confidence : bool, optional
            add confidence area to the regression line
        n_sigma : int, optional
            number of standard deviations for the confidence area.
        fs : int, optional
            fontsize

        """
        import matplotlib.pyplot as plt
        if ax is None:
            plt.figure()
            ax=plt.subplot(111)

        self.tree.root._v=0
        for n in self.tree.get_nonterminals(order='preorder'):
            for c in n:
                c._v = n._v + self.branch_value(c)

        tips = self.tree.get_terminals()
        internal = self.tree.get_nonterminals()

        # get values of terminals
        xi = np.array([self.tip_value(n) for n in tips])
        yi = np.array([n._v for n in tips])
        ind = np.array([n.bad_branch  if hasattr(n, 'bad_branch') else False for n in tips])
        if add_internal:
            xi_int = np.array([n.numdate for n in internal])
            yi_int = np.array([n._v for n in internal])
            ind_int = np.array([n.bad_branch  if hasattr(n, 'bad_branch') else False  for n in internal])

        if regression:
            # plot regression line
            t_mrca = -regression['intercept']/regression['slope']
            if add_internal:
                time_span = np.max(xi_int[~ind_int]) - np.min(xi_int[~ind_int])
                x_vals = np.array([max(np.min(xi_int[~ind_int]), t_mrca) - 0.1*time_span, np.max(xi[~ind])+0.05*time_span])
            else:
                time_span = np.max(xi[~ind]) - np.min(xi[~ind])
                x_vals = np.array([max(np.min(xi[~ind]), t_mrca) - 0.1*time_span, np.max(xi[~ind]+0.05*time_span)])

            # plot confidence interval
            if confidence and 'cov' in regression:
                x_vals = np.linspace(x_vals[0], x_vals[1], 100)
                y_vals = regression['slope']*x_vals + regression['intercept']
                dev = n_sigma*np.array([np.sqrt(regression['cov'][:2,:2].dot(np.array([x, 1])).dot(np.array([x,1]))) for x in x_vals])
                dev_slope = n_sigma*np.sqrt(regression['cov'][0,0])
                ax.fill_between(x_vals, y_vals-dev, y_vals+dev, alpha=0.2)
                dp = np.array([regression['intercept']/regression['slope']**2,-1./regression['slope']])
                dev_rtt = n_sigma*np.sqrt(regression['cov'][:2,:2].dot(dp).dot(dp))

            else:
                dev_rtt = None
                dev_slope = None

            ax.plot(x_vals, regression['slope']*x_vals + regression['intercept'],
                    label = r"$y=\alpha + \beta t$"+"\n"+
                            r"$\beta=$%1.2e"%(regression["slope"])
                            + ("+/- %1.e"%dev_slope if dev_slope else "") +
                            "\nroot date: %1.1f"%(-regression['intercept']/regression['slope']) +
                            ("+/- %1.2f"%dev_rtt if dev_rtt else ""))


        ax.scatter(xi[~ind], yi[~ind], label=("tips" if add_internal else None))
        if ind.sum():
            try:
                # note: this is treetime specific
                tmp_x = np.array([np.mean(n.raw_date_constraint) if n.raw_date_constraint else None
                                  for n in self.tree.get_terminals()])
                ax.scatter(tmp_x[ind], yi[ind], label="ignored tips", c='r')
            except:
                pass
        if add_internal:
            ax.scatter(xi_int[~ind_int], yi_int[~ind_int], label="internal nodes")

        ax.set_ylabel('root-to-tip distance', fontsize=fs)
        ax.set_xlabel('date', fontsize=fs)
        ax.ticklabel_format(useOffset=False)
        ax.tick_params(labelsize=fs*0.8)
        ax.set_ylim([0, 1.1*np.max(yi)])
        plt.tight_layout()
        plt.legend(fontsize=fs*0.8)
Plot root-to-tip distance vs time as a basic time-tree diagnostic

        Parameters
        ----------
        add_internal : bool, optional
            add internal nodes. this will only work if the tree has been dated already
        ax : None, optional
            an matplotlib axis to plot into. if non provided, a new figure is opened
        regression : None, optional
            a dict containing parameters of a root-to-tip vs time regression as
            returned by the function base_regression
        confidence : bool, optional
            add confidence area to the regression line
        n_sigma : int, optional
            number of standard deviations for the confidence area.
        fs : int, optional
            fontsize