Plot PCoA principal coordinates scaled by the relative abundances of
    otu_name.
491289
def plot_PCoA(cat_data, otu_name, unifrac, names, colors, xr, yr, outDir,
              save_as, plot_style):
    
    fig = plt.figure(figsize=(14, 8))
    ax = fig.add_subplot(111)

    for i, cat in enumerate(cat_data):
        plt.scatter(cat_data[cat]["pc1"], cat_data[cat]["pc2"], cat_data[cat]["size"],
                    color=colors[cat], alpha=0.85, marker="o", edgecolor="black",
                    label=cat)
    lgnd = plt.legend(loc="best", scatterpoints=3, fontsize=13)
    for i in range(len(colors.keys())):
        lgnd.legendHandles[i]._sizes = [80]  # Change the legend marker size manually
    plt.title(" ".join(otu_name.split("_")), style="italic")
    plt.ylabel("PC2 (Percent Explained Variance {:.3f}%)".format(float(unifrac["varexp"][1])))
    plt.xlabel("PC1 (Percent Explained Variance {:.3f}%)".format(float(unifrac["varexp"][0])))
    plt.xlim(round(xr[0]*1.5, 1), round(xr[1]*1.5, 1))
    plt.ylim(round(yr[0]*1.5, 1), round(yr[1]*1.5, 1))
    if plot_style:
        gu.ggplot2_style(ax)
        fc = "0.8"
    else:
        fc = "none"
    fig.savefig(os.path.join(outDir, "_".join(otu_name.split())) + "." + save_as,
                facecolor=fc, edgecolor="none", format=save_as,
                bbox_inches="tight", pad_inches=0.2)
    plt.close(fig)
830796
def create_scree_plot(data, o_filename, options):
    
    # Importing plt
    mpl.use("Agg")
    import matplotlib.pyplot as plt
    plt.ioff()

    # Computing the cumulative sum
    cumul_data = np.cumsum(data)

    # Creating the figure and axes
    fig, axes = plt.subplots(2, 1, figsize=(8, 16))

    # The title of the figure
    fig.suptitle(options.scree_plot_title, fontsize=16, weight="bold")

    fig.subplots_adjust(hspace=0.27, top=0.93, bottom=0.06)

    # Modifying the spines
    for axe in axes:
        axe.xaxis.set_ticks_position("bottom")
        axe.yaxis.set_ticks_position("left")
        axe.spines["top"].set_visible(False)
        axe.spines["right"].set_visible(False)
        axe.spines["left"].set_position(("outward", 9))
        axe.spines["bottom"].set_position(("outward", 9))

    # First, plotting the eigenvalues
    axes[0].set_title("Scree Plot", weight="bold")
    axes[0].set_xlabel("Component")
    axes[0].set_ylabel("Eigenvalue")
    axes[0].plot(np.arange(len(data)) + 1, data, marker="o", c="#0099CC",
                 mec="#0099CC", ls="-", lw=2, clip_on=False)

    # Then, plotting the annotation
    for i in range(len(data)):
        axes[0].annotate(np.round(data[i], 3),
                         xy=(np.arange(len(data))[i] + 1, data[i]),
                         xytext=(1, 10), textcoords="offset points",
                         ha="left", va="bottom",
                         bbox=dict(boxstyle="round,pad=0.5", fc="#FFFFFF"))

    # Next plot the cumulative values
    axes[1].set_title(("Cumulative explained variance "
                       "(max={:.3f})".format(np.sum(data))), weight="bold")
    axes[1].set_xlabel("Component")
    axes[1].set_ylabel("Cumulative explained variance")
    axes[1].axhline(np.sum(data) * 0.8, ls="--", lw="2", c="#999999")
    axes[1].plot(np.arange(len(data)) + 1, cumul_data, marker="o", c="#CC0000",
                 mec="#CC0000", mfc="#CC0000", ls="-", lw=2, clip_on=False)

    # Then, plotting the annotation
    for i in range(len(data)):
        axes[1].annotate(np.round(cumul_data[i], 3),
                         xy=(np.arange(len(data))[i] + 1, cumul_data[i]),
                         xytext=(1, -10), textcoords="offset points",
                         ha="left", va="top",
                         bbox=dict(boxstyle="round,pad=0.5", fc="#FFFFFF"))

    # Saving the file
    plt.savefig(o_filename, dpi=300)
    plt.close(fig)
354029
def _plot_posterior_op(
    values,
    var_name,
    selection,
    ax,
    bw,
    linewidth,
    bins,
    kind,
    point_estimate,
    round_to,
    credible_interval,
    ref_val,
    rope,
    ax_labelsize,
    xt_labelsize,
    **kwargs
):  # noqa: D202
    """"""

    def format_as_percent(x, round_to=0):
        return "{0:.{1:d}f}%".format(100 * x, round_to)

    def display_ref_val():
        if ref_val is None:
            return
        elif isinstance(ref_val, dict):
            val = None
            for sel in ref_val.get(var_name, []):
                if all(
                    k in selection and selection[k] == v for k, v in sel.items() if k != "ref_val"
                ):
                    val = sel["ref_val"]
                    break
            if val is None:
                return
        elif isinstance(ref_val, Number):
            val = ref_val
        else:
            raise ValueError(
                "Argument `ref_val` must be None, a constant, or a "
                'dictionary like {"var_name": {"ref_val": (lo, hi)}}'
            )
        less_than_ref_probability = (values < val).mean()
        greater_than_ref_probability = (values >= val).mean()
        ref_in_posterior = "{} <{:g}< {}".format(
            format_as_percent(less_than_ref_probability, 1),
            val,
            format_as_percent(greater_than_ref_probability, 1),
        )
        ax.axvline(val, ymin=0.05, ymax=0.75, color="C1", lw=linewidth, alpha=0.65)
        ax.text(
            values.mean(),
            plot_height * 0.6,
            ref_in_posterior,
            size=ax_labelsize,
            color="C1",
            weight="semibold",
            horizontalalignment="center",
        )

    def display_rope():
        if rope is None:
            return
        elif isinstance(rope, dict):
            vals = None
            for sel in rope.get(var_name, []):
                if all(k in selection and selection[k] == v for k, v in sel.items() if k != "rope"):
                    vals = sel["rope"]
                    break
            if vals is None:
                return
        elif len(rope) == 2:
            vals = rope
        else:
            raise ValueError(
                "Argument `rope` must be None, a dictionary like"
                '{"var_name": {"rope": (lo, hi)}}, or an'
                "iterable of length 2"
            )

        ax.plot(
            vals,
            (plot_height * 0.02, plot_height * 0.02),
            lw=linewidth * 5,
            color="C2",
            solid_capstyle="round",
            zorder=0,
            alpha=0.7,
        )
        text_props = {"size": ax_labelsize, "horizontalalignment": "center", "color": "C2"}
        ax.text(vals[0], plot_height * 0.2, vals[0], weight="semibold", **text_props)
        ax.text(vals[1], plot_height * 0.2, vals[1], weight="semibold", **text_props)

    def display_point_estimate():
        if not point_estimate:
            return
        if point_estimate not in ("mode", "mean", "median"):
            raise ValueError("Point Estimate should be in ('mode','mean','median')")
        if point_estimate == "mean":
            point_value = values.mean()
        elif point_estimate == "mode":
            if isinstance(values[0], float):
                density, lower, upper = _fast_kde(values, bw=bw)
                x = np.linspace(lower, upper, len(density))
                point_value = x[np.argmax(density)]
            else:
                point_value = mode(values.round(round_to))[0][0]
        elif point_estimate == "median":
            point_value = np.median(values)
        point_text = "{}={:.{}f}".format(point_estimate, point_value, round_to)

        ax.text(
            point_value,
            plot_height * 0.8,
            point_text,
            size=ax_labelsize,
            horizontalalignment="center",
        )

    def display_hpd():
        hpd_intervals = hpd(values, credible_interval=credible_interval)
        ax.plot(
            hpd_intervals,
            (plot_height * 0.02, plot_height * 0.02),
            lw=linewidth * 2,
            color="k",
            solid_capstyle="round",
        )
        ax.text(
            hpd_intervals[0],
            plot_height * 0.07,
            hpd_intervals[0].round(round_to),
            size=ax_labelsize,
            horizontalalignment="center",
        )
        ax.text(
            hpd_intervals[1],
            plot_height * 0.07,
            hpd_intervals[1].round(round_to),
            size=ax_labelsize,
            horizontalalignment="center",
        )
        ax.text(
            (hpd_intervals[0] + hpd_intervals[1]) / 2,
            plot_height * 0.3,
            format_as_percent(credible_interval) + " HPD",
            size=ax_labelsize,
            horizontalalignment="center",
        )

    def format_axes():
        ax.yaxis.set_ticks([])
        ax.spines["top"].set_visible(False)
        ax.spines["right"].set_visible(False)
        ax.spines["left"].set_visible(False)
        ax.spines["bottom"].set_visible(True)
        ax.xaxis.set_ticks_position("bottom")
        ax.tick_params(
            axis="x", direction="out", width=1, length=3, color="0.5", labelsize=xt_labelsize
        )
        ax.spines["bottom"].set_color("0.5")

    if kind == "kde" and values.dtype.kind == "f":
        plot_kde(
            values,
            bw=bw,
            fill_kwargs={"alpha": kwargs.pop("fill_alpha", 0)},
            plot_kwargs={"linewidth": linewidth},
            ax=ax,
            rug=False,
        )
    else:
        if bins is None:
            if values.dtype.kind == "i":
                xmin = values.min()
                xmax = values.max()
                bins = range(xmin, xmax + 2)
                ax.set_xlim(xmin - 0.5, xmax + 0.5)
            else:
                bins = "auto"
        kwargs.setdefault("align", "left")
        kwargs.setdefault("color", "C0")
        ax.hist(values, bins=bins, alpha=0.35, **kwargs)

    plot_height = ax.get_ylim()[1]

    format_axes()
    display_hpd()
    display_point_estimate()
    display_ref_val()
    display_rope()
578670
def plot_porosity_profile(self, fig=None):
        
        """
        if hasattr(self, '_fiber_image') is False:
            logger.warning('This method only works when a fiber image exists')
            return

        l = sp.asarray(sp.shape(self._fiber_image))
        px = sp.zeros(l[0])
        py = sp.zeros(l[1])
        pz = sp.zeros(l[2])

        for x in sp.arange(l[0]):
            px[x] = sp.sum(self._fiber_image[x, :, :])
            px[x] /= sp.size(self._fiber_image[x, :, :])
        for y in sp.arange(l[1]):
            py[y] = sp.sum(self._fiber_image[:, y, :])
            py[y] /= sp.size(self._fiber_image[:, y, :])
        for z in sp.arange(l[2]):
            pz[z] = sp.sum(self._fiber_image[:, :, z])
            pz[z] /= sp.size(self._fiber_image[:, :, z])

        if fig is None:
            fig = plt.figure()
        ax = fig.gca()
        plots = []
        plots.append(plt.plot(sp.arange(l[0])/l[0], px, 'r', label='x'))
        plots.append(plt.plot(sp.arange(l[1])/l[1], py, 'g', label='y'))
        plots.append(plt.plot(sp.arange(l[2])/l[2], pz, 'b', label='z'))
        plt.xlabel('Normalized Distance')
        plt.ylabel('Porosity')
        handles, labels = ax.get_legend_handles_labels()
        ax.legend(handles, labels, loc=1)
        plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
        return fig
75888
def plot_ossos_discoveries(ax, discoveries, plot_discoveries,
                           plot_colossos=False, split_plutinos=False):
    
    fc = ['b', '#E47833', 'k']
    alpha = [0.85, 0.6, 1.]
    marker = ['o', 'd']
    size = [7, 25]

    plottable = []   # Which blocks' discoveries to include?
    for d in discoveries:
        for n in plot_discoveries:
            if d['object'].startswith(n):  # can for sure be better, but this hack works. Need to get where going
                plottable.append(d)

    # Hack to get in the O15BD objects
    # directory_name = '/Users/bannisterm/Dropbox/OSSOS/measure3/ossin/D_tmp/'
    # kbos = parsers.ossos_discoveries(directory_name, all_objects=False, data_release=None)
    # for kbo in kbos:
    #     plottable_kbo = {'RAdeg': kbo.discovery.coordinate.ra.to_string(unit=units.degree, sep=':'),
    #                      'dist': kbo.orbit.distance.value}
    #     plottable.append(plottable_kbo)

    if plot_colossos:
        fainter = []
        colossos = []
        for n in plottable:
            if n['object'] in parameters.COLOSSOS:
                colossos.append(n)
            else:
                fainter.append(n)
        plot_ossos_points(fainter, ax, marker[0], size[0], fc[0], alpha[1], 1)
        plot_ossos_points(colossos, ax, marker[1], size[1], fc[2], alpha[2], 2)
    elif split_plutinos:
        # plutino_index = np.where((plottable['cl'] == 'res') & (plottable['j'] == 3) & (plottable['k'] == 2))
        raise NotImplementedError
    else:
        plot_ossos_points(plottable, ax, marker[0], size[0], fc[0], alpha[0], 2)

    return
802976
def composite_correlation(df, size=(12, 8)):
    
    fig = plt.figure(1, size)
    ax1 = plt.subplot2grid((2, 2), (0, 0))
    ax2 = plt.subplot2grid((2, 2), (0, 1))
    ax3 = plt.subplot2grid((2, 2), (1, 0))
    ax4 = plt.subplot2grid((2, 2), (1, 1))
    chemistry = ["V1", "V2", "V2.5", float("nan")]
    colors = sns.color_palette("Set2", 8)
    color_map = dict(zip(chemistry, colors))

    age_label = "Chronological age (yr)"
    ax1.scatter(df["hli_calc_age_sample_taken"], df["teloLength"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax1.set_ylim(0, 15)
    ax1.set_ylabel("Telomere length (Kb)")

    ax2.scatter(df["hli_calc_age_sample_taken"], df["ccn.chrX"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax2.set_ylim(1.8, 2.1)
    ax2.set_ylabel("ChrX copy number")

    ax4.scatter(df["hli_calc_age_sample_taken"], df["ccn.chrY"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax4.set_ylim(0.8, 1.1)
    ax4.set_ylabel("ChrY copy number")

    ax3.scatter(df["hli_calc_age_sample_taken"], df["TRA.PPM"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax3.set_ylim(0, 250)
    ax3.set_ylabel("$TCR-\\alpha$ deletions (count per million reads)")

    from matplotlib.lines import Line2D
    legend_elements = [Line2D([0], [0], marker='.', color='w', label=chem,
                          markerfacecolor=color, markersize=16) \
                        for (chem, color) in zip(chemistry, colors)[:3]]
    for ax in (ax1, ax2, ax3, ax4):
        ax.set_xlabel(age_label)
        ax.legend(handles=legend_elements, loc="upper right")

    plt.tight_layout()
    root = fig.add_axes((0, 0, 1, 1))
    labels = ((.02, .98, "A"),
              (.52, .98, "B"),
              (.02, .5, "C"),
              (.52, .5, "D"))
    panel_labels(root, labels)
    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()
979496
def _parse_coords(coord_lines):
        
        paras = {}
        var_pattern = re.compile(r"^([A-Za-z]+\S*)[\s=,]+([\d\-\.]+)$")
        for l in coord_lines:
            m = var_pattern.match(l.strip())
            if m:
                paras[m.group(1).strip("=")] = float(m.group(2))

        species = []
        coords = []
        # Stores whether a Zmatrix format is detected. Once a zmatrix format
        # is detected, it is assumed for the remaining of the parsing.
        zmode = False
        for l in coord_lines:
            l = l.strip()
            if not l:
                break
            if (not zmode) and GaussianInput._xyz_patt.match(l):
                m = GaussianInput._xyz_patt.match(l)
                species.append(m.group(1))
                toks = re.split(r"[,\s]+", l.strip())
                if len(toks) > 4:
                    coords.append([float(i) for i in toks[2:5]])
                else:
                    coords.append([float(i) for i in toks[1:4]])
            elif GaussianInput._zmat_patt.match(l):
                zmode = True
                toks = re.split(r"[,\s]+", l.strip())
                species.append(toks[0])
                toks.pop(0)
                if len(toks) == 0:
                    coords.append(np.array([0, 0, 0]))
                else:
                    nn = []
                    parameters = []
                    while len(toks) > 1:
                        ind = toks.pop(0)
                        data = toks.pop(0)
                        try:
                            nn.append(int(ind))
                        except ValueError:
                            nn.append(species.index(ind) + 1)
                        try:
                            val = float(data)
                            parameters.append(val)
                        except ValueError:
                            if data.startswith("-"):
                                parameters.append(-paras[data[1:]])
                            else:
                                parameters.append(paras[data])
                    if len(nn) == 1:
                        coords.append(np.array([0, 0, parameters[0]]))
                    elif len(nn) == 2:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        axis = [0, 1, 0]
                        op = SymmOp.from_origin_axis_angle(coords1, axis,
                                                           angle, False)
                        coord = op.operate(coords2)
                        vec = coord - coords1
                        coord = vec * bl / np.linalg.norm(vec) + coords1
                        coords.append(coord)
                    elif len(nn) == 3:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        coords3 = coords[nn[2] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        dih = parameters[2]
                        v1 = coords3 - coords2
                        v2 = coords1 - coords2
                        axis = np.cross(v1, v2)
                        op = SymmOp.from_origin_axis_angle(
                            coords1, axis, angle, False)
                        coord = op.operate(coords2)
                        v1 = coord - coords1
                        v2 = coords1 - coords2
                        v3 = np.cross(v1, v2)
                        adj = get_angle(v3, axis)
                        axis = coords1 - coords2
                        op = SymmOp.from_origin_axis_angle(
                            coords1, axis, dih - adj, False)
                        coord = op.operate(coord)
                        vec = coord - coords1
                        coord = vec * bl / np.linalg.norm(vec) + coords1
                        coords.append(coord)

        def _parse_species(sp_str):
            """
            The species specification can take many forms. E.g.,
            simple integers representing atomic numbers ("8"),
            actual species string ("C") or a labelled species ("C1").
            Sometimes, the species string is also not properly capitalized,
            e.g, ("c1"). This method should take care of these known formats.
            """
            try:
                return int(sp_str)
            except ValueError:
                sp = re.sub(r"\d", "", sp_str)
                return sp.capitalize()

        species = [_parse_species(sp) for sp in species]

        return Molecule(species, coords)
46547
def branchScale(self):
        """"""
        bs = -(self.prx * scipy.diagonal(self.Prxy, axis1=1, axis2=2)
                ).sum() * self.mu / float(self.nsites)
        assert bs > 0
        return bs
42414
def row_coordinates(self, X):
        """"""
        utils.validation.check_is_fitted(self, 'V_')

        _, row_names, _, _ = util.make_labels_and_names(X)

        if isinstance(X, pd.SparseDataFrame):
            X = X.to_coo().astype(float)
        elif isinstance(X, pd.DataFrame):
            X = X.to_numpy()

        if self.copy:
            X = X.copy()

        # Normalise the rows so that they sum up to 1
        if isinstance(X, np.ndarray):
            X = X / X.sum(axis=1)[:, None]
        else:
            X = X / X.sum(axis=1)

        return pd.DataFrame(
            data=X @ sparse.diags(self.col_masses_.to_numpy() ** -0.5) @ self.V_.T,
            index=row_names
        )
827462
def major_axis(points):
    
    U, S, V = np.linalg.svd(points)
    axis = util.unitize(np.dot(S, V))
    return axis
1008567
def aligned_residuals(pca):
    
    A = pca.rotated()
    fig, axes = P.subplots(2,1,
            sharex=True, frameon=False)
    fig.subplots_adjust(hspace=0, wspace=0.1)
    kw = dict(c="#555555", s=40, alpha=0.5)

    #lengths = attitude.pca.singular_values[::-1]
    lengths = (A[:,i].max()-A[:,i].min() for i in range(3))

    titles = (
        "Long cross-section (axis 3 vs. axis 1)",
        "Short cross-section (axis 3 vs. axis 2)")

    for title,ax,(a,b) in zip(titles,axes,
            [(0,2),(1,2)]):

        seaborn.regplot(A[:,a], A[:,b], ax=ax)
        ax.text(0,1,title,
            verticalalignment='top',
            transform=ax.transAxes)
        ax.autoscale(tight=True)
        for spine in ax.spines.itervalues():
            spine.set_visible(False)
    ax.set_xlabel("Meters")
    return fig
1234201
def plotPowers(ax, params, popkeys, dataset, linestyles, linewidths, transient=200, SCALING_POSTFIX='', markerstyles=None):
    '''
    '''
    
    colors = phlp.get_colors(len(popkeys))
    
    depth = params.electrodeParams['z']
    zpos = np.r_[params.layerBoundaries[:, 0],
                 params.layerBoundaries[-1, 1]]

    for i, layer in enumerate(popkeys):
        f = h5py.File(os.path.join(params.populations_path,
                                   '%s_population_%s' % (layer, dataset) + SCALING_POSTFIX + '.h5' ))
        ax.semilogx(f['data'].value[:, transient:].var(axis=1), depth,
                 color=colors[i],
                 ls=linestyles[i],
                 lw=linewidths[i],
                 marker=None if markerstyles is None else markerstyles[i],
                 markersize=2.5,
                 markerfacecolor=colors[i],
                 markeredgecolor=colors[i],
                 label=layer,
                 clip_on=True
                 )
    
        f.close()
    
    f = h5py.File(os.path.join(params.savefolder, '%ssum' % dataset + SCALING_POSTFIX + '.h5' ))
    ax.plot(f['data'].value[:, transient:].var(axis=1), depth,
                 'k', label='SUM', lw=1.25, clip_on=False)
    
    f.close()

    ax.set_yticks(zpos)
    ax.set_yticklabels([])
    #ax.set_xscale('log')
    try: # numticks arg only exists for latest matplotlib version
        ax.xaxis.set_major_locator(plt.LogLocator(base=10,
                                    subs=np.linspace(-10, 10, 2), numticks=6))
    except:
        ax.xaxis.set_major_locator(plt.LogLocator(base=10,
                                    subs=np.linspace(-10, 10, 2)))
    ax.xaxis.set_minor_locator(plt.LogLocator(base=10, subs=[1.]))
    ax.axis('tight')
175818
def fit(self, data):
        
        jmodel = callMLlibFunc("fitPCA", self.k, data)
        return PCAModel(jmodel)
981759
def gp_ccX():
    """"""
    inDir, outDir = getWorkDirs()
    data, alldata = OrderedDict(), None
    for infile in os.listdir(inDir):
        # get key and import data
        key = os.path.splitext(infile)[0].replace('_', '/')
        data_import = np.loadtxt(open(os.path.join(inDir, infile), 'rb'))
        # convert to log10 vs log10 plot, z=log10(y) => dz=0.434*dy/y
        data_import[:,3] = 0.434*(data_import[:,3]/data_import[:,1])
        data_import[:,(0,1)] = np.log10(data_import[:,(0,1)])
        data_import[:,2] = 0
        # fill dictionary
        data[key] = data_import
        alldata = data[key] if alldata is None else np.vstack((alldata, data[key]))
    # fit linear part first
    lindata = alldata[alldata[:,0]>2.5]
    m = (lindata[-1,1]-lindata[0,1])/(lindata[-1,0]-lindata[0,0])
    t = lindata[0,1] - m * lindata[0,0]
    popt1, pcov = curve_fit(
        linear, lindata[:,0], lindata[:,1], p0=[m, t],
        sigma=lindata[:,3], absolute_sigma=True
    )
    # fit full range
    popt2, pcov = curve_fit(
        lambda x, c, d: fitfunc(x, popt1[0], popt1[1], c, d),
        alldata[:,0], alldata[:,1], sigma=alldata[:,3], absolute_sigma=True,
    )
    popt = np.hstack((popt1, popt2))
    model = lambda x: fitfunc(x, *popt)
    # calculate mean standard deviation of data from parameterization
    yfit = np.array([model(x) for x in alldata[:,0]])
    stddev = 1.5*np.sqrt( # multiple of "sigma"
        np.average((alldata[:,1]-yfit)**2, weights=1./alldata[:,3])
    )
    print 'stddev = %.2g' % stddev
    errorband = np.array([[x, model(x), 0, 0, stddev] for x in np.linspace(1,4)])
    # make plot
    fitdata = np.array([[x, model(x), 0, 0, 0] for x in np.linspace(1,4)])
    par_names = ['a', 'b', 'c', 'd']
    energies = [19.6, 27, 39, 62.4, 200]
    labels = dict(
        ('%s = %.3g' % (par_name, popt[i]), [3.3, 3-i*0.2, True])
        for i,par_name in enumerate(par_names)
    )
    ccX_vals = [10**model(np.log10(energy)) for energy in energies]
    ccX = [' '.join([
        '%g GeV:' % energy,
        '({})'.format(ufloat(ccX_vals[i], stddev/0.434*ccX_vals[i])),
        '{/Symbol \155}b'
    ]) for i,energy in enumerate(energies)]
    print ccX
    #labels.update(dict(
    #    (cc, [1+i*0.5, 4.5+(i%2+1)*0.2, True]) for i,cc in enumerate(ccX)
    #))
    make_plot(
        data = [errorband] + data.values() + [fitdata],
        properties = [
            'with filledcurves lt 1 lw 5 pt 0 lc %s' % default_colors[8]
        ] + [
            'lc %s lw 4 lt 1 pt 18 ps 1.5' % (default_colors[i])
            for i in xrange(len(data))
        ] + ['with lines lc 0 lw 4 lt 1'],
        titles = [''] + data.keys() + ['y = ax+b - e^{-cx+d}'],
        xlabel = 'x = log_{10}[{/Symbol \326}s_{NN} (GeV)]',
        ylabel = 'y = log_{10}[{/Symbol \163}@_{c@^{/=18-}c}^{NN} ({/Symbol \155}b)]',
        name = os.path.join(outDir, 'ccX'),
        size = '11.4in,8.3in', xr = [1, 4], yr = [0.5,4.5],
        key = ['bottom right', 'nobox', 'width -5'],
        tmargin = 0.98, rmargin = 0.99, bmargin = 0.13,
        lines = dict(
            ('y=%f' % (np.log10(energy)), 'lc 0 lw 2 lt 2') for energy in energies
        ), labels = labels,
    )
578839
def plotConvergenceByObjectMultiColumn(results, objectRange, columnRange):
  
  ########################################################################
  #
  # Accumulate all the results per column in a convergence array.
  #
  # Convergence[c,o] = how long it took it to converge with f unique features
  # and c columns.

  convergence = numpy.zeros((max(columnRange), max(objectRange) + 1))
  for r in results:
    if r["numColumns"] in columnRange:
      convergence[r["numColumns"] - 1, r["numObjects"]] += r["convergencePoint"]

  convergence /= numTrials

  # print "Average convergence array=", convergence

  ########################################################################
  #
  # Create the plot. x-axis=
  plt.figure()
  plotPath = os.path.join("plots", "convergence_by_object_multicolumn.jpg")

  # Plot each curve
  legendList = []
  colorList = ['r', 'b', 'g', 'm', 'c', 'k', 'y']

  for i in range(len(columnRange)):
    c = columnRange[i]
    print "columns={} objectRange={} convergence={}".format(
      c, objectRange, convergence[c-1,objectRange])
    if c == 1:
      legendList.append('1 column')
    else:
      legendList.append('{} columns'.format(c))
    plt.plot(objectRange, convergence[c-1,objectRange],
             color=colorList[i])

  # format
  plt.legend(legendList, loc="upper left", prop={'size':10})
  plt.xlabel("Number of objects in training set")
  plt.xticks(range(0,max(objectRange)+1,10))
  plt.yticks(range(0,int(convergence.max())+2))
  plt.ylabel("Average number of touches")
  plt.title("Object recognition with multiple columns (unique features = 5)")

    # save
  plt.savefig(plotPath)
  plt.close()
844636
def slamdunkPCAPlot (self):
        

        pconfig = {
            'id': 'slamdunk_pca',
            'title': 'Slamdunk: PCA',
            'xlab': 'PC1',
            'ylab': 'PC2',
            'tt_label': 'PC1 {point.x:.2f}: PC2 {point.y:.2f}'
        }

        self.add_section (
            name = 'PCA (T&gt;C based)',
            anchor = 'slamdunk_PCA',
            description = """This plot shows the principal components of samples based
                        on the distribution of reads with T&gt;C conversions within UTRs
                        (see the <a href="http://t-neumann.github.io/slamdunk/docs.html#summary" target="_blank">slamdunk docs</a>).""",
            plot = scatter.plot(self.PCA_data, pconfig)
        )
953136
def compare_variant_type_plot(data):
    
    keys = OrderedDict()
    keys['snps'] = {'name': 'SNPs', 'color': '#7cb5ec'}
    keys['indels'] = {'name': 'InDels', 'color': '#90ed7d'}
    keys['multiallelic_snps'] = {'name': 'multi-allelic SNP', 'color': 'orange'}
    keys['complex_indels'] = {'name': 'Complex InDels', 'color': '#8085e9'}

    total_variants = dict()
    known_variants = dict()
    novel_variants = dict()
    for s_name, values in data.items():
        total_variants[s_name] = {
            'snps': values['TOTAL_SNPS'],
            'indels': values['TOTAL_INDELS'],
            'multiallelic_snps': values['TOTAL_MULTIALLELIC_SNPS'],
            'complex_indels': values['TOTAL_COMPLEX_INDELS'],
        }

        known_variants[s_name] = {
            'snps': values['NUM_IN_DB_SNP'],
            'indels': int(values['TOTAL_INDELS'])-int(values['NOVEL_INDELS']),
            'multiallelic_snps': values['NUM_IN_DB_SNP_MULTIALLELIC'],
            'complex_indels': values['NUM_IN_DB_SNP_COMPLEX_INDELS'],
        }

        novel_variants[s_name] = {
            'snps': values['NOVEL_SNPS'],
            'indels': values['NOVEL_INDELS'],
            'multiallelic_snps': int(values['TOTAL_MULTIALLELIC_SNPS'])-int(values['NUM_IN_DB_SNP_MULTIALLELIC']),
            'complex_indels': int(values['TOTAL_COMPLEX_INDELS'])-int(values['NUM_IN_DB_SNP_COMPLEX_INDELS']),
        }

    plot_conf = {
        'id': 'picard_variantCallingMetrics_variant_type',
        'title': 'Picard: Variants Called',
        'ylab': 'Counts of Variants',
        'hide_zero_cats': False,
        'data_labels': [
            {'name': 'Total'},
            {'name': 'Known'},
            {'name': 'Novel'}
        ],
    }
    return bargraph.plot(data=[total_variants, known_variants, novel_variants],
                         cats=[keys, keys, keys],
                         pconfig=plot_conf)