def plot_results(t_hr, I_mA, avg_t_hr, lifetime_hr, lifetime_error = None,
                 show_plot = False):
    """"""
    
    # Set show_plot to True to stop the code here for examining the plot
    if show_plot:
        fig = plt.figure()
            
        ax1 = fig.add_subplot(111)
        h1 = ax1.plot(t_hr, I_mA, '.-', label='Data', linewidth=0.5)
        ax1.set_xlabel('Time [hr]')
        ax1.set_ylabel('Beam Current [mA]')
            
        ax2 = ax1.twinx()
        
        if lifetime_error is not None :
            h2 = ax2.errorbar(avg_t_hr, lifetime_hr,
                              yerr = lifetime_error,
                              fmt = '-', ecolor = 'r',
                              color = 'r', linestyle = '-', marker = '.', 
                              label = 'Lifetime')
        else :
            h2 = ax2.plot(avg_t_hr, lifetime_hr,
                          color = 'r', linestyle = '-', marker = '.',
                          label = 'Lifetime')
            
        ax2.set_ylabel('Lifetime [hr]')
            
        plt.figlegend( (h1[0], h2[0]), 
                       (h1[0].get_label(), h2[0].get_label()), 
                       loc='upper right')

        plt.show()  

  def plot(self,yl='',yr='',x='s',idx=slice(None),clist='k r b g c m',lattice=True,newfig=True):
    self._clist=clist.split()
    if newfig:
      f=_p.figure()
    else:
      f=_p.gcf()
    sp=_p.subplot(111)
    _p.subplots_adjust(right=0.72)
    _p.setp( sp.yaxis, visible=False)
    xd=getattr(self,x)[idx]
    out=qdplot(figure=f,subplot=sp,plots=[],lines=[],legends=[],xaxis=xd)
    if lattice:
      self._lattice(out,idx,['kn0l','angle'],"#a0ffa0",'Bend h')
      self._lattice(out,idx,['ks0l'],"#ffa0a0",'Bend v')
      self._lattice(out,idx,['kn1l','k1l'],"#a0a0ff",'Quad')
      self._lattice(out,idx,['hkick'],"#e0a0e0",'Kick h')
      self._lattice(out,idx,['vkick'],"#a0e0e0",'Kick v')
    for i in yl.split():
      self._column(out,idx,i,'left')
    for i in yr.split():
      self._column(out,idx,i,'right')

    _p.xlabel(_mylbl(axlabel,x))

    _p.xlim(min(xd),max(xd))
    _p.figlegend(out.lines,out.legends,'upper right')
    _p.grid()
    _p.draw()
    del self._clist
    if hasattr(self,'_colAxleft'): delattr(self,'_colAxleft')
    if hasattr(self,'_colAxright'): delattr(self,'_colAxright')
    self.currplot=out
    return out

def plot(x, y, yci, methods, title, xlabel, ylabel, filename):
  fig = pylab.figure()

  ax = pylab.gca()
  for method in methods:
    print method, y(method), yci(method)
    lines = ax.errorbar(x(method), y(method), yci(method), fmt=markers[method], label=legends[method], mfc='none', markersize=15, capsize=5)
  
  # easiest way to plot brute force
  """
  method = 'brute'
  lines = ax.errorbar(range(1,4), [0,] * 3, [0,] * 3, fmt=markers[method], label=legends[method], mfc='none', markersize=15, capsize=5)
  """

  ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f'))

  pylab.title(title)
  pylab.xlabel(xlabel)
  pylab.ylabel(ylabel)
  #pylab.ylim([0, 0.2])
  pylab.gcf().subplots_adjust(bottom=0.2, left=0.2)
  fig.savefig(filename + ".pdf", dpi=300, format="pdf")
  
  figLegend = pylab.figure(figsize = (6, 4))
  pylab.figlegend(*ax.get_legend_handles_labels(), loc = 'upper left')
  figLegend.savefig("legend.pdf", dpi=300, format="pdf")
  
  pylab.close()

  def update(self):
    if callable(self.pre):
      self.pre()
    _p.ioff()
    self.lines=[]
    self.legends=[]
#    self.figure.lines=[]
#    self.figure.patches=[]
#    self.figure.texts=[]
#    self.figure.images = []
    self.figure.legends = []

    if self.lattice:
      self.lattice.patches=[]
      self._lattice(['k0l','kn0l','angle'],"#a0ffa0",'Bend h')
      self._lattice(['ks0l'],"#ffa0a0",'Bend v')
      self._lattice(['kn1l','k1l'],"#a0a0ff",'Quad')
      self._lattice(['hkick'],"#e0a0e0",'Kick h')
      self._lattice(['vkick'],"#a0e0e0",'Kick v')
    if self.left:
      self.left.lines=[]
      for i in self.yl:
        self._column(i,self.left,self.color[i])
    if self.right:
      self.right.lines=[]
      for i in self.yr:
        self._column(i,self.right,self.color[i])
    _p.xlabel(_mylbl(axlabel,self.x))
    self.figure.gca().set_xlim(min(self.xaxis[self.idx]),max(self.xaxis[self.idx]))
    _p.figlegend(self.lines,self.legends,'upper right')
    _p.grid(True)
#    self.figure.canvas.mpl_connect('button_release_event',self.button_press)
#    self.figure.canvas.mpl_connect('pick_event',self.pick)
    _p.ion()
    _p.draw()

def plotMatches(imgsources, refsources, matches, wcs, W, H, prefix,
                saveplot=True, format='png'):
    plt.clf()

    # Image sources
    ix = np.array([s.getXAstrom() for s in imgsources])
    iy = np.array([s.getYAstrom() for s in imgsources])
    iflux = np.array([s.getPsfFlux() for s in imgsources])
    I = np.argsort(-iflux)
    # First 200: red dots
    II = I[:200]
    p1 = plt.plot(ix[II], iy[II], 'r.', zorder=10)
    # Rest: tiny dots
    II = I[200:]
    p2 = plt.plot(ix[II], iy[II], 'r.', markersize=1, zorder=9)

    # Ref sources:
    # Only getRa() (not getRaAstrom(), getRaObject()) is non-zero.

    rx,ry = [],[]
    for r in refsources:
        xy = wcs.skyToPixel(r.getRaDec())
        rx.append(xy[0])
        ry.append(xy[1])
    rx = np.array(rx)
    ry = np.array(ry)
    p3 = plt.plot(rx, ry, 'bo', mec='b', mfc='none', markersize=6, zorder=20)

    x,y = [],[]
    dx,dy = [],[]
    for m in matches:
        x0,x1 = m.first.getXAstrom(), m.second.getXAstrom()
        y0,y1 = m.first.getYAstrom(), m.second.getYAstrom()
        #plt.plot([x0, x1], [y0, y1], 'g.-')
        x.append(x0)
        y.append(y0)
        dx.append(x1-x0)
        dy.append(y1-y0)
    #plt.plot(x, y, 's', mec='g', mfc='none', markersize=5)
    p4 = plt.plot(x, y, 'o', mec='g', mfc='g', alpha=0.5, markersize=8, zorder=5)
    p5 = plt.quiver(x, y, dx, dy, angles='xy', scale=30., zorder=30)
    plt.axis('scaled')
    plt.axis([0, W, 0, H])
    #print p1, p2, p3, p4, p5

    plt.figlegend((p1, p2, p3, p4), #, p5),
              ('Image sources (brightest 200)',
               'Image sources (rest)',
               'Reference sources',
               'Matches',),
              'center right',
              numpoints=1,
              prop=FontProperties(size='small'))

    fn = prefix + '-matches.' + format
    return _output(fn, format, saveplot)

def plot(cluster, filename="plot.png", func=lambda a: a.id, 
         plot_title='', cmap='Paired', filter=lambda a: True, 
         draw_legend=False, radius='2.25', sym=None):
    ac = rc.load(cluster)
    clusters = rc.load_clusters(cluster)

    p = get_population()
    pops = [0 for i in range(30000)]
    cluster_pops = [] 
    cluster_pop_max = []
    for clust in range(len(clusters)):
        cluster_pops.append(pops[:])
        cluster_pop_max.append([0,0,-1,0])

    for clust,agents in enumerate(clusters):
        for agent in agents:
            a = Agent(agent)
            for i in range(a.birth, a.death):
                cluster_pops[clust][i] += 1
                if cluster_pop_max[clust][2] == -1:
                    cluster_pop_max[clust][2] = i
                if i > cluster_pop_max[clust][3]:
                    cluster_pop_max[clust][3] = i
                if cluster_pops[clust][i] > cluster_pop_max[clust][0]:
                    cluster_pop_max[clust][0] = cluster_pops[clust][i]
                    cluster_pop_max[clust][1] = i
    
    lines=[]
    for i,clust in enumerate(cluster_pops):
        lines.append(pylab.plot(range(30000),clust, 
                         label=("%d: k=%d" % (i, len(clusters[i]))),
                         color=pylab.cm.Paired(float(i)/len(clusters))))

    if draw_legend:
        pylab.figlegend(lines, ["%d: k=%d" % (i, len(clust))
                                    for i,clust in enumerate(clusters)], 
                        'center right', 
                        ncol=((len(clusters)/35)+1), prop=dict(size=6))
    else:
        print "should not draw!!!"

    title = r"Cluster Population ($\epsilon$ = %s, %d clusters)" % (radius, len(clusters))
    pylab.title(title, weight='black')
    pylab.xlabel("Time", weight='bold')
    pylab.ylabel("Population Size", weight='bold')
    if sym is not None:
        pylab.figtext(0,.954, '(%s)' % sym, size=6, weight='black')

    pylab.savefig(filename, dpi=300)

    print 'cluster, totalPop, start, peak, stop, maxPop'
    for clust,agents in enumerate(clusters):
        print clust, len(agents), cluster_pop_max[clust][2], cluster_pop_max[clust][1], cluster_pop_max[clust][3]+1, cluster_pop_max[clust][0]

 def piechart(self, domain_type):
     pie_file = self.output[0]
     pylab.figure(figsize=(5,5))
     
     probs,labels = ([],[])
     for type in DOMAIN_TYPES:
         mean_coverage = mean(domain_type[type]) if len(domain_type[type])>0 else 0.0
         debug(self.family_name,type,mean_coverage)
         probs.append(mean_coverage)
         labels.append(type)#+":%0.2f" % mean_coverage)
     
     #ax = pylab.axes([0.6, 0.6, 0.4, 0.4])
     explode = [0.05 for i in xrange(len(DOMAIN_TYPES))]
     patches, texts = pylab.pie(probs,explode=None,labels=None,shadow=False,colors=DOMAIN_COLORS)
     pylab.figlegend(patches, labels, "lower left", fancybox=True, markerscale=0.2)
     pylab.title(self.family_name)
     pylab.savefig(pie_file)

def plot_categories(categories, category_names=None, distribution_names=None, numrows=None, numcols=None):
    "Plot the categories as sub-plots."

    import pylab as P
    assert 2 == len(categories.shape)

    # determine number of rows and columns
    num_dists = categories.shape[0]
    num_categories = categories.shape[1]
    if not numrows and not numcols:
        numrows = int(numpy.sqrt(num_categories))
    if not numrows:
        numrows = num_categories / numcols
        if num_categories % numcols:
            numrows += 1
    elif not numcols:
        numcols = num_categories / numrows
        if num_categories % numrows:
            numcols += 1

    # plot each category
    max, min = categories.max(), categories.min()
    ind = numpy.arange(num_dists)    # the x locations for the groups
    # the width of the bars: can also be len(x) sequence
    width = 1.
    for c in range(num_categories):
        # overlaps, subplot(111) is killed
        splt = P.subplot(numrows, numcols, c + 1)
        handles = [
            P.bar((i,), (x,), width, color=color)
            for color, i, x
            in zip(pylab_utils.simple_colours, ind, categories[:, c])
        ]
        splt.set_ylim(min, max)
        if category_names:
            P.title(category_names[c])
        if c % numcols:
            splt.set_yticklabels(())
        if distribution_names:
            P.xticks(ind + width / 2., distribution_names)
        else:
            splt.set_xticklabels(())
        splt.set_xticks(tuple())
    if distribution_names:
        P.figlegend(handles, distribution_names, 'upper right')

def plot(cluster_file, filename="plot.png", func=lambda a: a.id, 
         plot_title='', cmap='Paired', filter=lambda a: True, 
         draw_legend=False, radius='2.25', sym=None, run_dir='../run/'):
    """
    Creates a line plot showing population per cluster
    """

    # retrieve cluster data from cluster file
    clusters = rc.load_clusters(cluster_file)

    # grab cluster population
    p = get_population(run_dir=run_dir)

    lines=[]
    for cluster, agents in enumerate(clusters):
        pop_by_time = list(repeat(0, 30000))
        for agent in agents:
            a = Agent(agent)
            for i in range(a.birth, a.death):
                pop_by_time[i] += 1
    
        lines.append(
            pylab.plot(range(30000), pop_by_time, 
                       label=("%d: k=%d" % (i, len(clusters[cluster]))),
                       color=pylab.cm.Paired(float(cluster)/len(clusters))))

    if draw_legend:
        pylab.figlegend(lines, ["%d: k=%d" % (i, len(cluster))
                                    for i,cluster in enumerate(clusters)], 
                        'center right', 
                        ncol=((len(clusters)/35)+1), prop=dict(size=6))

    title = r"Cluster Population ($\epsilon$ = %s, %d clusters)" % (radius, len(clusters))
    pylab.title(title, weight='black')
    pylab.xlabel("Time", weight='bold')
    pylab.ylabel("Population Size", weight='bold')
    if sym is not None:
        pylab.figtext(0,.954, '(%s)' % sym, size=6, weight='black')

    pylab.savefig(filename, dpi=300)

def do_plot(data, legend, fn_png):
    pylab.clf()
    pylab.figure(0, (6,8))
    #pylab.axes([left, bottom, width, height])
    # in fractional coordinates
    pylab.axes([0.20, 0.070, 0.76, 0.84])
    handles = []
    labels = []
    for key, symbol, label in legend:
        x = []
        y = []
        for index, lot in enumerate(lots_list):
            av = data.get((key, lot.key))
            if av is not None:
                x.append(av)
                y.append(-index)
        handle = pylab.plot(x, y, symbol)
        handles.append(handle)
        labels.append(label)
    pos = -numpy.arange(len(lots_list))
    ylabels = [lot.key for lot in lots_list]
    pylab.yticks(pos, ylabels, size="small")
    pylab.xticks(size="small")
    #pylab.xlabel(r"Geometric mean of $k_{\theory}/k_{\exp}$", size="small")
    pylab.xlabel(r"$\exp(RMSD(\ln(k_{\theory}) - \ln($k_{\exp})))$", size="small")
    pylab.semilogx()
    pylab.axvline(0.1, color="k", zorder=-5)
    pylab.axvline(1.0, color="k", zorder=-5)
    pylab.axvline(10.0, color="k", zorder=-5)
    for y in 0, -5, -10, -15, -20, -25, -30, -35, -40:
        pylab.axhline(y, color="k", alpha=0.2, lw=10, zorder=-5)
    if len(legend) > 1:
        legend = pylab.figlegend(handles, labels, (0.22,0.915), numpoints=1, handlelength=0, prop={"size":"small"})
        frame = legend.get_frame()
        frame.set_lw(0.0)
    #pylab.xlim(1e-8, 1e5)
    pylab.xlim(1e0, 1e5)
    pylab.ylim(pos[-1]-0.5, 0.5)
    pylab.savefig(fn_png)

    pylab.semilogy()
    pylab.fill(
        numpy.concatenate((invtemps, invtemps[::-1])),
        numpy.concatenate((experimental_k*10, experimental_k[::-1]/10)),
        "k", alpha=0.2, lw=0,
    )
    pylab.xticks(
        1/numpy.array([300, 350, 400, 450, 500, 550, 600], float),
        ["300", "350", "400", "450", "500", "550", "600"],
    )
    pylab.title(title)
    pylab.xlabel("T [K]")
    pylab.ylabel("k [(m**3/mol)/s]")
    pylab.ylim(1e-8,1e7)
    legend = pylab.figlegend(
        lines, labels, (0.07,0.06), ncol=3, prop={"size":11},
        handlelength=3, labelspacing=0.1, columnspacing=1
    )
    #legend.get_frame().set_linewidth(0)
    legend.get_frame().set_fill(True)
    legend.get_frame().set_alpha(0.5)
    pylab.savefig(fn_img % "rates")

    pylab.clf()
    lines = []
    labels = []
    line = pylab.plot([experimental_Ea], [experimental_A], color="k", marker="o", ms=11, mew=2, lw=0, ls=" ")
    lines.append(line)
    labels.append("experiment")
    for lot in lots_list:
        if lot.spin == "ROS":
            label = "ro" + lot.label

#.........这里部分代码省略.........
                    yerr = np.array([yerr for a in x])

                plots = []
                labels = []
                for i, name in enumerate(y.rnames):
                    if aggregate != None:
                        yerr = yerrs[i].flatten()
                    
                    labels.append(name[0][1])

                    if all([_isfloat(a) for a in x]):
                        plots.append(
                            axs[-1].errorbar(x, y[i].flatten(), yerr)[0])
                        
                        if xmin == 'AUTO' and xmax == 'AUTO':
                            xmin , xmax = axs[-1].get_xlim()
                            xran = xmax - xmin
                            xmin = xmin - .05*xran
                            xmax = xmax + .05*xran
                            
                        axs[-1].plot([xmin, xmax], [0.,0.], 'k:')
                        
                    else : # categorical x axis
                        plots.append(
                            axs[-1].errorbar(
                                _xrange(len(x)), y[i].flatten(), yerr)[0])
                        
                        pylab.xticks(_xrange(len(x)), x)
                        xmin = - 0.5
                        xmax = len(x) - 0.5
                        
                        axs[-1].plot([xmin, xmax], [0., 0.], 'k:')

                pylab.figlegend(plots, labels, loc=1,
                                labelsep=.005,
                                handlelen=.01,
                                handletextsep=.005)

            test['y'].append(y)
            if yerr == None:
                test['yerr'].append([])
            else:
                test['yerr'].append(yerr)
            test['xmins'].append(xmin)
            test['xmaxs'].append(xmax)

            #  8.2 Add subplot title
            ######################################################
            if rlevels == [1] and clevels == [1]:
                title = ''
                
            elif rlevels == [1]:
                title = _str(clevel)
                
            elif clevels == [1]:
                title = _str(rlevel)
                
            else:
                title = '%s = %s, %s = %s' \
                        % (sepyplots, _str(rlevel),
                           sepxplots, _str(rlevel))
                
            pylab.title(title, fontsize='medium')
            test['subplot_titles'].append(title)

            #  8.3 Format the subplot

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('log_dir', nargs='*')
    parser.add_argument('--n_cols', type=int, default=4)
    parser.add_argument('--smoothing', type=float, default=0.6)
    args = parser.parse_args()

    if not args.log_dir:
        raise Exception("Specify at least one log directory")
    for d in args.log_dir:
        if not os.path.exists(d):
            raise Exception(f"Directory {d} does not exist")

    print("Reading events...")
    events_by_log_dir = read_all_events(args.log_dir)

    i = len(os.path.commonprefix(args.log_dir))
    len_longest_suffix = len(os.path.commonprefix([s[::-1] for s in args.log_dir]))
    if len_longest_suffix > 0:
        j = -len_longest_suffix
    else:
        j = None

    print("Plotting...")
    lines, labels = [], []
    plot_n = 1
    axes = {}
    n_graphs = len(set([key for events in events_by_log_dir.values() for key in events.keys()])) + 1
    subplot_dims = (int(ceil(n_graphs / args.n_cols)), args.n_cols)
    figure(figsize=(6 * subplot_dims[1], 4 * subplot_dims[0]))
    for log_dir_n, (log_dir, events) in enumerate(events_by_log_dir.items()):
        print(log_dir)
        color = f"C{log_dir_n % 10}"
        label = log_dir[i:j]
        for key, events in events.items():
            if key not in axes:
                axes[key] = subplot(*(subplot_dims + (plot_n,)))
                plot_n += 1
            ax = axes[key]
            timestamps, values = list(zip(*events))
            relative_timestamps = [t - timestamps[0] for t in timestamps]
            relative_timestamps_hours = [t / 3600 for t in relative_timestamps]
            line = plot_values(ax, values, relative_timestamps_hours, 'Value', 'Time (hours)', key,
                               label=label, color=color, smoothing=args.smoothing)
            lines.append(line)
            labels.append(label)

    unique_lines = []
    unique_labels = []
    for line, label in zip(lines, labels):
        if label not in unique_labels:
            unique_labels.append(label)
            unique_lines.append(line)
    figlegend(unique_lines, unique_labels, loc='upper center')
    tight_layout(rect=[0, 0, 1, 0.90])

    print("Saving figure...")
    # Removing trailing slashes
    normed_log_dirs = [os.path.normpath(d) for d in args.log_dir]
    if len(args.log_dir) == 1:
        out_filename = os.path.basename(normed_log_dirs[0])
    else:
        longest_common_prefix = os.path.commonprefix(normed_log_dirs)
        longest_common_suffix = os.path.commonprefix([d[::-1] for d in normed_log_dirs])[::-1]
        out_filename = longest_common_prefix + '*' + longest_common_suffix + '.png'
    savefig(out_filename)

        P.close()

    # do AUC bar-chart
    #P.rcParams['xtick.direction'] = 'out'
    def x(a, b, m):
        return b*(len(methods)+1) + m
    P.figure(figsize=(14,6))
    for a, auc in enumerate(('AUC', 'AUC50')):
        ax = P.subplot(1, 2, a+1)
        xlocs = []
        xlabels = []
        for b, bg in enumerate(harness.options.backgrounds):
            rects = P.bar(
                [x(a, b, m) for m in xrange(len(methods))],
                [aucs[bg][auc][method] for method in methods],
                color=[colors[method] for method in methods],
                width=1.
            )
            xlocs.append(x(a, b, len(methods)/2.))
            xlabels.append(bg)
        P.xticks(xlocs, xlabels)
        P.ylim(0,1)
        P.title(auc)
        for tl in ax.get_xticklines():
            tl.set_visible(False)
    #P.text(.5, .9, '%s' % fragment, horizontalalignment='center')
    P.figlegend(rects, methods, loc='upper right')
    P.savefig(os.path.join(options.results_dir, 'AUC-%s.png' % fragment))
    P.savefig(os.path.join(options.results_dir, 'AUC-%s.eps' % fragment))
    P.close()

#.........这里部分代码省略.........


    # Add the box and whiskers
    galaxies2 = filter(lambda galaxy: galaxy.ston_J > 30., galaxies)
    galaxies2 = pyl.asarray(galaxies2)
    x = [galaxy.Mass for galaxy in galaxies2]
    ll = 8.5
    ul= 12
    #bins_x =pyl.linspace(ll, ul, 7)
    #bins_x =pyl.arange(8.5, 12.5, 0.5)
    bins_x =pyl.array([8.5, 9., 9.5, 10., 10.5, 11., 12.])
    grid = []

    for i in range(bins_x.size-1):
        xmin = bins_x[i]
        xmax = bins_x[i+1]
        cond=[cond1 and cond2 for cond1, cond2 in zip(x>=xmin, x<xmax)]

        grid.append(galaxies2.compress(cond))
    icd = []
    for i in range(len(grid)):
        icd.append([galaxy.ICD_JH*100 for galaxy in grid[i]])
    #bp2 = f1s2.boxplot(icd, positions=pyl.delete(bins_x,-1)+0.25, sym='')
    pbp(f1s3, icd, indexer=list(pyl.delete(bins_x,-1)+0.25))
    pbp(f1s4, icd, indexer=list(pyl.delete(bins_x,-1)+0.25))

    # Finish Plot
    # Tweak colors on the boxplot
    #pyl.setp(bp1['boxes'], lw=2)
    #pyl.setp(bp1['whiskers'], lw=2)
    #pyl.setp(bp1['medians'], lw=2)
    #pyl.setp(bp2['boxes'], lw=2)
    #pyl.setp(bp2['whiskers'], lw=2)
    #pyl.setp(bp2['medians'], lw=2)
    #pyl.setp(bp['fliers'], color='#8CFF6F', marker='+')

    #f1s1.axvspan(7.477, 9, facecolor='#FFFDD0', ec='None', zorder=0)
    #f1s1.axvspan(11, 12, facecolor='#FFFDD0', ec='None', zorder=0)
    #f1s2.axvspan(7.477, 9, facecolor='#FFFDD0', ec='None', zorder=0)
    #f1s2.axvspan(11, 12, facecolor='#FFFDD0', ec='None', zorder=0)

    f1s1.set_xlim(8,12)
    f1s2.set_xlim(8,12)
    f1s3.set_xlim(8,12)
    f1s4.set_xlim(8,12)

    f1s1.set_ylim(-10,50)
    f1s2.set_ylim(0,15)
    f1s3.set_ylim(-5,12)
    f1s4.set_ylim(-1,3)

    f1s1.set_xticks([8,9,10,11,12])
    f1s1.set_xticklabels([])
    f1s2.set_xticks([8,9,10,11,12])
    f1s2.set_xticklabels([])
    f1s3.set_xticks([8,9,10,11,12])
    f1s4.set_xticks([8,9,10,11,12])

    f1s4.set_yticks([-1, 0, 1, 2, 3])

    f1s3.set_xlabel(r"Log Mass ($M_{\odot})$")
    f1s1.set_ylabel(r"$\xi[i_{775},H_{160}]$ (%)")
    f1s4.set_xlabel(r"Log Mass ($M_{\odot})$")
    f1s3.set_ylabel(r"$\xi[J_{125},H_{160}]$ (%)")

    import matplotlib.font_manager
    line1 = pyl.Line2D([], [], marker='o', mfc='0.8', mec='0.8', markersize=8,
        linewidth=0)
    line2 = pyl.Line2D([], [], marker='s', mec='blue', mfc='None',
        markersize=10, linewidth=0, markeredgewidth=2)
    line3 = pyl.Line2D([], [], color='r', linewidth=2)
    prop = matplotlib.font_manager.FontProperties(size='small')
    pyl.figlegend((line1, line2, line3), ('Data', 'Quartiles',
        'Medians'), 'lower center', prop=prop, ncol=3)

    from matplotlib.patches import ConnectionPatch
    xy = (12, 15)
    xy2 = (8, 15)
    con = ConnectionPatch(xyA=xy, xyB=xy2, coordsA='data', coordsB='data',
        axesA=f1s1, axesB=f1s2)
    xy = (12, 0)
    xy2 = (8, 0)
    con2 = ConnectionPatch(xyA=xy, xyB=xy2, coordsA='data', coordsB='data',
        axesA=f1s1, axesB=f1s2)
    f1s1.add_artist(con)
    f1s1.add_artist(con2)

    xy = (12, 3)
    xy2 = (8, 3)
    con = ConnectionPatch(xyA=xy, xyB=xy2, coordsA='data', coordsB='data',
        axesA=f1s3, axesB=f1s4)
    xy = (12, -1)
    xy2 = (8, -1)
    con2 = ConnectionPatch(xyA=xy, xyB=xy2, coordsA='data', coordsB='data',
        axesA=f1s3, axesB=f1s4)
    f1s3.add_artist(con)
    f1s3.add_artist(con2)

    pyl.draw()
    pyl.show()

def plot_icd_vs_mass():
    galaxies = pickle.load(open('galaxies.pickle','rb'))

    arrow_up = [[0,0], [-1,-1], [0,0], [0,-2], [0,0], [1,-1], [0,0]]

    # Make figure
    f1 = pyl.figure(1, figsize=(4,6))
    f1s1 = f1.add_subplot(211)
    f1s2 = f1.add_subplot(212)

    for galaxy in galaxies:
        if galaxy.ston_I > 30. and galaxy.ICD_IH != None:
            f1s1.scatter(galaxy.z, galaxy.ICD_IH * 100, c='0.8',
                marker='o', s=25, edgecolor='0.8')

        if galaxy.ston_V > 30. and galaxy.ICD_VJ != None:
            f1s2.scatter(galaxy.z, galaxy.ICD_VJ * 100, c='0.8',
                marker='o', s=25, edgecolor='0.8')

    # Add the box and whiskers
    galaxies2 = filter(lambda galaxy: galaxy.ston_I > 30., galaxies)
    galaxies2 = pyl.asarray(galaxies2)
    x = [galaxy.z for galaxy in galaxies2]
    ll = 1.5
    ul= 4
    bins_x =pyl.arange(ll, ul, 0.5)
    grid = []

    for i in range(bins_x.size-1):
        xmin = bins_x[i]
        xmax = bins_x[i+1]
        cond=[cond1 and cond2 for cond1, cond2 in zip(x>=xmin, x<xmax)]

        grid.append(galaxies2.compress(cond))
    icd = []
    for i in range(len(grid)):
        icd.append([galaxy.ICD_IH*100 for galaxy in grid[i]])
    from boxplot_percentile import percentile_box_plot as pbp
    pbp(f1s1, icd, indexer=list(pyl.delete(bins_x,-1)+0.25), whisker_top=None,
            whisker_bottom=None)

    # Add the box and whiskers
    galaxies = pickle.load(open('galaxies.pickle','rb'))
    galaxies2 = filter(lambda galaxy: galaxy.ston_V > 30., galaxies)
    galaxies2 = pyl.asarray(galaxies2)
    x = [galaxy.z for galaxy in galaxies2]
    ll = 1.5
    ul= 4
    bins_x =pyl.arange(ll, ul, 0.5)
    grid = []

    for i in range(bins_x.size-1):
        xmin = bins_x[i]
        xmax = bins_x[i+1]
        cond=[cond1 and cond2 for cond1, cond2 in zip(x>=xmin, x<xmax)]

        grid.append(galaxies2.compress(cond))
    icd = []
    for i in range(len(grid)):
        icd.append([galaxy.ICD_VJ*100 for galaxy in grid[i]])
    pbp(f1s2, icd, indexer=list(pyl.delete(bins_x,-1)+0.25), whisker_top=None,
            whisker_bottom=None)

    f1s1.set_xlim(1.5, 3.5)
    f1s2.set_xlim(1.5, 3.5)

    f1s1.set_ylim(-1, 10)
    f1s2.set_ylim(-1, 10)

    f1s1.set_xticklabels([])
    f1s2.set_xticks([1.5,2,2.5,3,3.5])

    f1s1.set_ylabel(r"$\xi[i_{775},H_{160}]$ (%)")
    f1s2.set_xlabel("Redshift")
    f1s2.set_ylabel(r"$\xi[V_{606},J_{125}]$ (%)")

    import matplotlib.font_manager
    line1 = pyl.Line2D([], [], marker='o', mfc='0.8', mec='0.8', markersize=8,
        linewidth=0)
    line2 = pyl.Line2D([], [], marker='s', mec='#348ABD', mfc='None',
        markersize=10, linewidth=0, markeredgewidth=2)
    line3 = pyl.Line2D([], [], color='#A60628', linewidth=2)
    prop = matplotlib.font_manager.FontProperties(size='small')
    pyl.figlegend((line1, line2, line3), ('Data', 'Quartiles',
        'Medians'), 'lower center', prop=prop, ncol=3)

    pyl.show()

#.........这里部分代码省略.........

            D_results_build[algorithm][i] = (t1 - t0)
            D_results_query[algorithm][i] = (t2 - t1)

    #------------------------------------------------------------
    # varying k
    k_results_build = dict([(alg, np.zeros(len(krange)))
                            for alg in algorithms])
    k_results_query = dict([(alg, np.zeros(len(krange)))
                            for alg in algorithms])

    X = get_data(N, DD, dataset)

    for i, kk in enumerate(krange):
        print("k = %i (%i out of %i)" % (kk, i + 1, len(krange)))
        for algorithm in algorithms:
            nbrs = neighbors.NearestNeighbors(n_neighbors=kk,
                                              algorithm=algorithm,
                                              leaf_size=leaf_size)
            t0 = time()
            nbrs.fit(X)
            t1 = time()
            nbrs.kneighbors(X)
            t2 = time()

            k_results_build[algorithm][i] = (t1 - t0)
            k_results_query[algorithm][i] = (t2 - t1)

    pl.figure(figsize=(8, 11))

    for (sbplt, vals, quantity,
         build_time, query_time) in [(311, Nrange, 'N',
                                      N_results_build,
                                      N_results_query),
                                     (312, Drange, 'D',
                                      D_results_build,
                                      D_results_query),
                                     (313, krange, 'k',
                                      k_results_build,
                                      k_results_query)]:
        ax = pl.subplot(sbplt, yscale='log')
        pl.grid(True)

        tick_vals = []
        tick_labels = []

        bottom = 10 ** np.min([min(np.floor(np.log10(build_time[alg])))
                               for alg in algorithms])

        for i, alg in enumerate(algorithms):
            xvals = 0.1 + i * (1 + len(vals)) + np.arange(len(vals))
            width = 0.8

            c_bar = pl.bar(xvals, build_time[alg] - bottom,
                           width, bottom, color='r')
            q_bar = pl.bar(xvals, query_time[alg],
                           width, build_time[alg], color='b')

            tick_vals += list(xvals + 0.5 * width)
            tick_labels += ['%i' % val for val in vals]

            pl.text((i + 0.02) / len(algorithms), 0.98, alg,
                    transform=ax.transAxes,
                    ha='left',
                    va='top',
                    bbox=dict(facecolor='w', edgecolor='w', alpha=0.5))

            pl.ylabel('Time (s)')

        ax.xaxis.set_major_locator(ticker.FixedLocator(tick_vals))
        ax.xaxis.set_major_formatter(ticker.FixedFormatter(tick_labels))

        for label in ax.get_xticklabels():
            label.set_rotation(-90)
            label.set_fontsize(10)

        title_string = 'Varying %s' % quantity

        descr_string = ''

        for s in 'NDk':
            if s == quantity:
                pass
            else:
                descr_string += '%s = %i, ' % (s, fiducial_values[s])

        descr_string = descr_string[:-2]

        pl.text(1.01, 0.5, title_string,
                transform=ax.transAxes, rotation=-90,
                ha='left', va='center', fontsize=20)

        pl.text(0.99, 0.5, descr_string,
                transform=ax.transAxes, rotation=-90,
                ha='right', va='center')

        pl.gcf().suptitle("%s data set" % dataset.capitalize(), fontsize=16)

    pl.figlegend((c_bar, q_bar), ('construction', 'N-point query'),
                 'upper right')

#.........这里部分代码省略.........
        pylab.ylabel(options.ytitle)

    lines = []
    # use dummy_lines to workaround a bug in errorbars that
    # causes the line styles to be set incorrectly.
    dummy_lines = []
    new_legend = []

    if options.error_bars:
        if options.error_bars == "interleaved":
            step_size = 2
            max_size = len(legend)
        elif options.error_bars == "blocked":
            step_size = 1
            max_size = (len(legend) - 1) / 2
    else:
        step_size = 1
        max_size = len(legend)

    if options.stacked:
        dummy_lines, new_legend = doStackedPlot(data, legend)
    else:
        nplotted = 0
        nskipped = 0
        for x in range(1, max_size, step_size):

            s = options.symbols[nplotted % len(options.symbols)]

            yvals = data[:, x]

            xvals = numpy.ma.masked_array(data[:, 0], numpy.ma.getmask(yvals))

            xvals = xvals.compressed()
            yvals = yvals.compressed()

            if len(xvals) == 0:
                E.warn("skipped empty column %i: %s" % (x, legend[x]))

            if options.error_bars == "interleaved":
                yerr = data[:, x + 1]
                yerr = yerr.compressed()
            else:
                yerr = None

            lines.append(pylab.errorbar(xvals,
                                        yvals,
                                        yerr=yerr,
                                        fmt=s))

            dummy_lines.append(pylab.plot(xvals,
                                          yvals,
                                          s))

            new_legend.append(legend[x])

            nplotted += 1

        E.info("nplotted=%i, nskipped=%i" % (nplotted, nskipped))

    if len(lines) == 0:
        E.Stop()
        return

    if options.legend_location != "none":
        pylab.figlegend(dummy_lines,
                        new_legend,
                        options.legend_location)

    if options.logscale:
        if "x" in options.logscale:
            pylab.gca().set_xscale('log')
        if "y" in options.logscale:
            pylab.gca().set_yscale('log')

    if options.xrange:
        pylab.xlim(options.xrange)

    if options.yrange:
        pylab.ylim(options.yrange)

    if options.function:
        xstart, xend = pylab.gca().get_xlim()
        increment = (xend - xstart) / 100.0
        exec("f = lambda x: %s" % options.function) in locals()
        xvals, yvals = [], []
        for x in range(0, 100):
            xvals.append(xstart)
            yvals.append(f(xstart))
            xstart += increment
        xvals.append(xstart)
        yvals.append(f(xstart))

        pylab.plot(xvals, yvals)

    if options.hardcopy:
        pylab.savefig(os.path.expanduser(options.hardcopy), dpi=options.dpi)
    else:
        pylab.show()

    E.Stop()

    sp = subplot(211) # switch to the first subplot
    cla() # clear the subplot
    title("Training Set") # set the subplot's title
    sp.set_autoscale_on( True ) # enable autoscaling
    targetline = plot(trnSequenceTarget,"r-") # plot the targets
    sp.set_autoscale_on( False ) # disable autoscaling
    outputline = plot(trnSequenceOutput,"b-") # plot the actual output


    # plot test data
    sp = subplot(212)
    cla()
    title("Test Set")
    sp.set_autoscale_on( True )
    plot(tstSequenceTarget,"r-")
    sp.set_autoscale_on( False )
    plot(tstSequenceOutput,"b-")

    # create a legend
    figlegend((targetline, outputline),('target','output'),('upper right'))

    # draw everything
    draw()


show()

    # calculate the nets output for train and the test data
    trnSequenceOutput = net.extrapolate(trnSequenceWashout, len(trnSequenceTarget))
    tstSequenceOutput = net.extrapolate(tstSequenceWashout, len(tstSequenceTarget))

    # plot training data
    sp = subplot(211)  # switch to the first subplot
    cla()  # clear the subplot
    title("Training Set")  # set the subplot's title
    sp.set_autoscale_on(True)  # enable autoscaling
    targetline = plot(trnSequenceTarget, "r-")  # plot the targets
    sp.set_autoscale_on(False)  # disable autoscaling
    outputline = plot(trnSequenceOutput, "b-")  # plot the actual output

    # plot test data
    sp = subplot(212)
    cla()
    title("Test Set")
    sp.set_autoscale_on(True)