def visualize(segmentation, expression, visualize=None, store=None, title=None, legend=False):
    notes = []
    onsets = []
    values = []
    param = ['Dynamics', 'Articulation', 'Tempo']
    converter = NoteList()
    converter.bpm = 100
    if not visualize:
        visualize = selectSubset(param)
    for segment, expr in zip(segmentation, expression):
        for note in segment:
            onsets.append(converter.ticks_to_milliseconds(note.on)/1000.0)
            values.append([expr[i] for i in visualize])
    import matplotlib.pyplot as plt
    fig = plt.figure(figsize=(12, 4))
    for i in visualize:
        plt.plot(onsets, [v[i] for v in values], label=param[i])
    plt.ylabel('Deviation')
    plt.xlabel('Score time (seconds)')
    if legend:
        plt.legend(bbox_to_anchor=(0., 1), loc=2, borderaxespad=0.)

    if title:
        plt.title(title)
    #dplot = fig.add_subplot(111)
    #sodplot = fig.add_subplot(111)
    #dplot.plot([i for i in range(len(deltas[0]))], deltas[0])
    #sodplot.plot([i for i in range(len(sodeltas[0]))], sodeltas[0])
    if store:
        fig.savefig('plots/{0}.png'.format(store))
    else:
        plt.show()

def showCumulOverlap(mode, modes, *args, **kwargs):
    """Show cumulative overlap using :func:`~matplotlib.pyplot.plot`.

    :type mode: :class:`.Mode`, :class:`.Vector`
    :arg modes: multiple modes
    :type modes: :class:`.ModeSet`, :class:`.ANM`, :class:`.GNM`, :class:`.PCA`
    """

    import matplotlib.pyplot as plt
    if not isinstance(mode, (Mode, Vector)):
        raise TypeError('mode must be NMA, ModeSet, Mode or Vector, not {0}'
                        .format(type(mode)))
    if not isinstance(modes, (NMA, ModeSet)):
        raise TypeError('modes must be NMA, ModeSet, or Mode, not {0}'
                        .format(type(modes)))
    cumov = (calcOverlap(mode, modes) ** 2).cumsum() ** 0.5
    if isinstance(modes, NMA):
        arange = np.arange(0.5, len(modes)+0.5)
    else:
        arange = modes.getIndices() + 0.5
    show = plt.plot(arange, cumov, *args, **kwargs)
    plt.title('Cumulative overlap with {0}'.format(str(mode)))
    plt.xlabel('{0} mode index'.format(modes))
    plt.ylabel('Cumulative overlap')
    plt.axis((arange[0]-0.5, arange[-1]+0.5, 0, 1))
    if SETTINGS['auto_show']:
        showFigure()
    return show

def showNormDistFunct(model, coords, *args, **kwargs):
    """Show normalized distance fluctuation matrix using 
    :func:`~matplotlib.pyplot.imshow`. By default, *origin=lower* 
    keyword  arguments are passed to this function, 
    but user can overwrite these parameters."""

    import math
    import matplotlib
    import matplotlib.pyplot as plt
    normdistfunct = model.getNormDistFluct(coords)

    if not 'origin' in kwargs:
        kwargs['origin'] = 'lower'
        
    matplotlib.rcParams['font.size'] = '14'
    fig = plt.figure(num=None, figsize=(10,8), dpi=100, facecolor='w')
    show = plt.imshow(normdistfunct, *args, **kwargs), plt.colorbar()
    plt.clim(math.floor(np.min(normdistfunct[np.nonzero(normdistfunct)])), \
                                           round(np.amax(normdistfunct),1))
    plt.title('Normalized Distance Fluctution Matrix')
    plt.xlabel('Indices', fontsize='16')
    plt.ylabel('Indices', fontsize='16')
    if SETTINGS['auto_show']:
        showFigure()
    return show

def showOverlapTable(modes_x, modes_y, **kwargs):
    """Show overlap table using :func:`~matplotlib.pyplot.pcolor`.  *modes_x*
    and *modes_y* are sets of normal modes, and correspond to x and y axes of
    the plot.  Note that mode indices are incremented by 1.  List of modes
    is assumed to contain a set of contiguous modes from the same model.

    Default arguments for :func:`~matplotlib.pyplot.pcolor`:

      * ``cmap=plt.cm.jet``
      * ``norm=plt.normalize(0, 1)``"""

    import matplotlib.pyplot as plt

    overlap = abs(calcOverlap(modes_y, modes_x))
    if overlap.ndim == 0:
        overlap = np.array([[overlap]])
    elif overlap.ndim == 1:
        overlap = overlap.reshape((modes_y.numModes(), modes_x.numModes()))

    cmap = kwargs.pop('cmap', plt.cm.jet)
    norm = kwargs.pop('norm', plt.normalize(0, 1))
    show = (plt.pcolor(overlap, cmap=cmap, norm=norm, **kwargs),
            plt.colorbar())
    x_range = np.arange(1, modes_x.numModes() + 1)
    plt.xticks(x_range-0.5, x_range)
    plt.xlabel(str(modes_x))
    y_range = np.arange(1, modes_y.numModes() + 1)
    plt.yticks(y_range-0.5, y_range)
    plt.ylabel(str(modes_y))
    plt.axis([0, modes_x.numModes(), 0, modes_y.numModes()])
    if SETTINGS['auto_show']:
        showFigure()
    return show

def showNormedSqFlucts(modes, *args, **kwargs):
    """Show normalized square fluctuations via :func:`~matplotlib.pyplot.plot`.
    """

    import matplotlib.pyplot as plt
    sqf = calcSqFlucts(modes)
    args = list(args)
    modesarg = []
    i = 0
    while i < len(args):
        if isinstance(args[i], (VectorBase, ModeSet, NMA)):
            modesarg.append(args.pop(i))
        else:
            i += 1
    show = [plt.plot(sqf/(sqf**2).sum()**0.5, *args,
                     label='{0}'.format(str(modes)), **kwargs)]
    plt.xlabel('Indices')
    plt.ylabel('Square fluctuations')
    for modes in modesarg:
        sqf = calcSqFlucts(modes)
        show.append(plt.plot(sqf/(sqf**2).sum()**0.5, *args,
                    label='{0}'.format(str(modes)), **kwargs))
    if SETTINGS['auto_show']:
        showFigure()
    return show

def scree_plot(pca_obj, fname=None): 
    '''
    Scree plot for variance & cumulative variance by component from PCA. 

    Arguments: 
        - pca_obj: a fitted sklearn PCA instance
        - fname: path to write plot to file

    Output: 
        - scree plot 
    '''   
    components = pca_obj.n_components_ 
    variance = pca.explained_variance_ratio_
    plt.figure()
    plt.plot(np.arange(1, components + 1), np.cumsum(variance), label='Cumulative Variance')
    plt.plot(np.arange(1, components + 1), variance, label='Variance')
    plt.xlim([0.8, components]); plt.ylim([0.0, 1.01])
    plt.xlabel('No. Components', labelpad=11); plt.ylabel('Variance Explained', labelpad=11)
    plt.legend(loc='best') 
    plt.tight_layout() 
    if fname is not None:
        plt.savefig(fname)
        plt.close() 
    else:
        plt.show() 
    return 

def draw_ranges_for_parameters(data, title='', save_path='./pictures/'):
  parameters = data.columns.values.tolist()

  # remove flight name parameter
  for idx, parameter in enumerate(parameters):
    if parameter == 'flight_name':
      del parameters[idx]

  flight_names = np.unique(data['flight_name'])

  print len(flight_names)

  for parameter in parameters:
    plt.figure()

    axis = plt.gca()

    # ax.set_xticks(numpy.arange(0,1,0.1))
    axis.set_yticks(flight_names)
    axis.tick_params(labelright=True)
    axis.set_ylim([94., 130.])
    plt.grid()

    plt.title(title)
    plt.xlabel(parameter)
    plt.ylabel('flight name')

    colors = iter(cm.rainbow(np.linspace(0, 1,len(flight_names))))

    for flight in flight_names:
      temp = data[data.flight_name == flight][parameter]

      plt.plot([np.min(temp), np.max(temp)], [flight, flight], c=next(colors), linewidth=2.0)
    plt.savefig(save_path+title+'_'+parameter+'.jpg')
    plt.close()

def plotTestData(tree):
	plt.figure()
	plt.axis([0,1,0,1])
	plt.xlabel("X axis")
	plt.ylabel("Y axis")
	plt.title("Green: Class1, Red: Class2, Blue: Class3, Yellow: Class4")
	for value in class1:
		plt.plot(value[0],value[1],'go')
	plt.hold(True)
	for value in class2:
		plt.plot(value[0],value[1],'ro')
	plt.hold(True)
	for value in class3:
		plt.plot(value[0],value[1],'bo')
	plt.hold(True)
	for value in class4:
		plt.plot(value[0],value[1],'yo')
	plotRegion(tree)
	for value in classPlot1:
		plt.plot(value[0],value[1],'g.',ms=3.0)
	plt.hold(True)
	for value in classPlot2:
		plt.plot(value[0],value[1],'r.', ms=3.0)
	plt.hold(True)
	for value in classPlot3:
		plt.plot(value[0],value[1],'b.', ms=3.0)
	plt.hold(True)
	for value in classPlot4:
		plt.plot(value[0],value[1],'y.', ms=3.0)
	plt.grid(True)
	plt.show()

def plotJ(J_history,num_iters):
    x = np.arange(1,num_iters+1)
    plt.plot(x,J_history)
    plt.xlabel(u"迭代次数",fontproperties=font) # 注意指定字体，要不然出现乱码问题
    plt.ylabel(u"代价值",fontproperties=font)
    plt.title(u"代价随迭代次数的变化",fontproperties=font)
    plt.show()

def plotISVar():
    plt.figure()
    plt.title('Variance minimization problem (call).\nVertical lines mark the minima.')
    for K in [0.6, 0.8, 1.0, 1.2]:
        theta = np.linspace(-0.6, 2)
        var = [BS.exactCallVar(K*s0, theta) for theta in theta]
        minth = theta[np.argmin(var)]
        line, = plt.plot(theta, var, label=str(K))
        plt.axvline(minth, color=line.get_color())

    plt.xlabel(r'$\theta$')
    plt.ylabel('call variance')
    plt.legend(title=r'$K/s_0$', loc='upper left')
    plt.autoscale(tight=True)

    plt.figure()
    plt.title('Variance minimization problem (put).\nVertical lines mark the minima.')
    for K in [0.8, 1.0, 1.2, 1.4]:
        theta = np.linspace(-2, 0.5)
        var = [BS.exactPutVar(K*s0, theta) for theta in theta]
        minth = theta[np.argmin(var)]
        line, = plt.plot(theta, var, label=str(K))
        plt.axvline(minth, color=line.get_color())

    plt.xlabel(r'$\theta$')
    plt.ylabel('put variance')
    plt.legend(title=r'$K/s_0$', loc='upper left')
    plt.autoscale(tight=True)

def plot_mpl_fig(): 
    rootdir = '/Users/catherinefielder/Documents/Research_Halos/HaloDetail'
    cs = []
    pops = []
    for subdir, dirs, files in os.walk(rootdir):
        head,tail = os.path.split(subdir)
        haloname = tail
        for file in files:
            if file.endswith('_columnsadded'):
                values = ascii.read(os.path.join(subdir, file), format = 'commented_header') #Get full path and access file
                host_c = values[1]['host_c']  
                cs = np.append(cs, host_c)                     
                pop = len(values['mvir(10)'])
                pops = np.append(pops, pop)
                print pop
                plt.loglog(host_c, pop, alpha=0.8,label = haloname)
        print "%s done. On to the next." %haloname
    #plt.xscale('log')
    #plt.yscale('log')
    plt.xlabel('Host Concentration')
    plt.ylabel('Nsat')
    plt.title('Abundance vs. Host Concentration', ha='center')
    #plt.legend(loc='best')
    spearman = scipy.stats.spearmanr(cs, pops)
    print spearman

def display(spectrum):
	template = np.ones(len(spectrum))

	#Get the plot ready and label the axes
	pyp.plot(spectrum)
	max_range = int(math.ceil(np.amax(spectrum) / standard_deviation))
	for i in range(0, max_range):
		pyp.plot(template * (mean + i * standard_deviation))
	pyp.xlabel('Units?')
	pyp.ylabel('Amps Squared')    
	pyp.title('Mean Normalized Power Spectrum')
	if 'V' in Options:
		pyp.show()
	if 'v' in Options:
		tokens = sys.argv[-1].split('.')
		filename = tokens[0] + ".png"
		input = ''
		if os.path.isfile(filename):
			input = input("Error: Plot file already exists! Overwrite? (y/n)\n")
			while input != 'y' and input != 'n':
				input = input("Please enter either \'y\' or \'n\'.\n")
			if input == 'y':
				pyp.savefig(filename) 
			else:
				print("Plot not written.")
		else:
			pyp.savefig(filename) 

def scatter(x, y, equal=False, xlabel=None, ylabel=None, xinvert=False, yinvert=False):
    """
    Plot a scatter with simple formatting options
    """
    plt.scatter(x, y, 200, color=[0.3, 0.3, 0.3], edgecolors="white", linewidth=1, zorder=2)
    sns.despine()
    if xlabel:
        plt.xlabel(xlabel)
    if ylabel:
        plt.ylabel(ylabel)
    if equal:
        plt.axes().set_aspect("equal")
        plt.plot([0, max([x.max(), y.max()])], [0, max([x.max(), y.max()])], color=[0.6, 0.6, 0.6], zorder=1)
        bmin = min([x.min(), y.min()])
        bmax = max([x.max(), y.max()])
        rng = abs(bmax - bmin)
        plt.xlim([bmin - rng * 0.05, bmax + rng * 0.05])
        plt.ylim([bmin - rng * 0.05, bmax + rng * 0.05])
    else:
        xrng = abs(x.max() - x.min())
        yrng = abs(y.max() - y.min())
        plt.xlim([x.min() - xrng * 0.05, x.max() + xrng * 0.05])
        plt.ylim([y.min() - yrng * 0.05, y.max() + yrng * 0.05])
    if xinvert:
        plt.gca().invert_xaxis()
    if yinvert:
        plt.gca().invert_yaxis()

def tuning(x, y, err=None, smooth=None, ylabel=None, pal=None):
    """
    Plot a tuning curve
    """
    if smooth is not None:
        xs, ys = smoothfit(x, y, smooth)
        plt.plot(xs, ys, linewidth=4, color="black", zorder=1)
    else:
        ys = asarray([0])
    if pal is None:
        pal = sns.color_palette("husl", n_colors=len(x) + 6)
        pal = pal[2 : 2 + len(x)][::-1]
    plt.scatter(x, y, s=300, linewidth=0, color=pal, zorder=2)
    if err is not None:
        plt.errorbar(x, y, yerr=err, linestyle="None", ecolor="black", zorder=1)
    plt.xlabel("Wall distance (mm)")
    plt.ylabel(ylabel)
    plt.xlim([-2.5, 32.5])
    errTmp = err
    errTmp[isnan(err)] = 0
    rng = max([nanmax(ys), nanmax(y + errTmp)])
    plt.ylim([0 - rng * 0.1, rng + rng * 0.1])
    plt.yticks(linspace(0, rng, 3))
    plt.xticks(range(0, 40, 10))
    sns.despine()
    return rng

def plotAlphas(datasetNames, sampleSizes, foldsSet, cvScalings, sampleMethods, fileNameSuffix): 
    """
    Plot the variation in the error with alpha for penalisation. 
    """
    for i, datasetName in enumerate(datasetNames): 
        #plt.figure(i)    
        
        
        for k in range(len(sampleMethods)):
            outfileName = outputDir + datasetName + sampleMethods[k] + fileNameSuffix + ".npz"
            data = numpy.load(outfileName)
    
            errors = data["arr_0"]
            meanMeasures = numpy.mean(errors, 0)
            
            foldInd = 4 
    
            for i in range(sampleSizes.shape[0]):
                plt.plot(cvScalings, meanMeasures[i, foldInd, 2:8], next(linecycler), label="m="+str(sampleSizes[i]))
                    
            plt.xlabel("Alpha")
            plt.ylabel('Error')
            xmin, xmax = cvScalings[0], cvScalings[-1]
            plt.xlim((xmin,xmax))

        
            plt.legend(loc="upper left")
    plt.show()

def plotIterationResult(train_err_list):
    x = range(1,len(train_err_list) + 1)
    fig = plt.figure()
    plt.plot(x,train_err_list)
    plt.xlabel('iterations')
    plt.ylabel('training error')
    plt.show()

def do_plot(mode, content, wide):
	global style
	style.apply(mode, content, wide)

	data = np.load("data/prr_AsAu_%s%s.npz"%(content, wide))

	AU, TAU = np.meshgrid(-data["Au_range_dB"], data["tau_range"])
	Zu = data["PRR_U"]
	Zs = data["PRR_S"]

	assert TAU.shape == AU.shape == Zu.shape, "The inputs TAU, AU, PRR_U must have the same shape for plotting!"

	plt.clf()

	if mode in ("sync",):
		# Plot the inverse power ratio, sync signal is stronger for positive ratios
		CSf = plt.contourf(TAU, AU, Zs, levels=(0.0, 0.2, 0.4, 0.6, 0.8, 0.9, 1.0), colors=("1.0", "0.75", "0.5", "0.25", "0.15", "0.0"), origin="lower")
		CS2 = plt.contour(CSf, colors = ("r",)*5+("w",), linewidths=(0.75,)*5+(1.0,), origin="lower", hold="on")
	else:
		CSf  = plt.contourf(TAU, AU, Zs, levels=(0.0, 0.2, 0.4, 0.6, 0.8, 0.9, 1.0), colors=("1.0", "0.75", "0.5", "0.25", "0.15", "0.0"), origin="lower")
		CS2f = plt.contour(CSf, levels=(0.0, 0.2, 0.4, 0.6, 0.8, 1.0), colors=4*("r",)+("w",), linewidths=(0.75,)*4+(1.0,), origin="lower", hold="on")
		#CS2f = plt.contour(TAU, -AU, Zu, levels=(0.9, 1.0), colors=("0.0",), linewidths=(1.0,), origin="lower", hold="on")
		if content in ("unif",):
			CSu  = plt.contourf(TAU, AU, Zu, levels=(0.2, 1.0), hatches=("////",), colors=("0.75",), origin="lower")
			CS2  = plt.contour(CSu, levels=(0.2,), colors = ("r",), linewidths=(1.0,), origin="lower", hold="on")

	style.annotate(mode, content, wide)

	plt.axis([data["tau_range"][0], data["tau_range"][-1], -data["Au_range_dB"][-1], -data["Au_range_dB"][0]])

	plt.ylabel(r"Signal power ratio ($\mathrm{SIR}$)", labelpad=2)
	plt.xlabel(r"Time offset $\tau$ ($/T$)", labelpad=2)

	plt.savefig("pdf/prrc2_%s_%s%s_z.pdf"%(mode, content, wide))

def plotErrorBars(dict_to_plot, x_lim, y_lim, xlabel, y_label, title, out_file, margin=[0.05, 0.05], loc=2):

    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(y_label)

    if y_lim is None:
        y_lim = [1 * float("Inf"), -1 * float("Inf")]

    max_val_seen_y = y_lim[1] - margin[1]
    min_val_seen_y = y_lim[0] + margin[1]
    print min_val_seen_y, max_val_seen_y
    max_val_seen_x = x_lim[1] - margin[0]
    min_val_seen_x = x_lim[0] + margin[0]
    handles = []
    for k in dict_to_plot:
        means, stds, x_vals = dict_to_plot[k]

        min_val_seen_y = min(min(np.array(means) - np.array(stds)), min_val_seen_y)
        max_val_seen_y = max(max(np.array(means) + np.array(stds)), max_val_seen_y)

        min_val_seen_x = min(min(x_vals), min_val_seen_x)
        max_val_seen_x = max(max(x_vals), max_val_seen_x)

        handle = plt.errorbar(x_vals, means, yerr=stds)
        handles.append(handle)
        print max_val_seen_y
    plt.xlim([min_val_seen_x - margin[0], max_val_seen_x + margin[0]])
    plt.ylim([min_val_seen_y - margin[1], max_val_seen_y + margin[1]])
    plt.legend(handles, dict_to_plot.keys(), loc=loc)
    plt.savefig(out_file)

    def plot(self, standardized=False, **kwargs):
        """
        standardized: standardize each estimated coefficient and confidence interval endpoints by the standard error of the estimate.

        """
        from matplotlib import pyplot as plt

        ax = kwargs.get('ax', None) or plt.figure().add_subplot(111)
        yaxis_locations = range(len(self.hazards_.columns))

        summary = self.summary
        lower_bound = self.confidence_intervals_.loc['lower-bound'].copy()
        upper_bound = self.confidence_intervals_.loc['upper-bound'].copy()
        hazards = self.hazards_.values[0].copy()

        if standardized:
            se = summary['se(coef)']
            lower_bound /= se
            upper_bound /= se
            hazards /= se

        order = np.argsort(hazards)
        ax.scatter(upper_bound.values[order], yaxis_locations, marker='|', c='k')
        ax.scatter(lower_bound.values[order], yaxis_locations, marker='|', c='k')
        ax.scatter(hazards[order], yaxis_locations, marker='o', c='k')
        ax.hlines(yaxis_locations, lower_bound.values[order], upper_bound.values[order], color='k', lw=1)

        tick_labels = [c + significance_code(p).strip() for (c, p) in summary['p'][order].iteritems()]
        plt.yticks(yaxis_locations, tick_labels)
        plt.xlabel("standardized coef" if standardized else "coef")
        return ax

  def build_hist(self, coverage, show=False, save=False, save_fn="max_hist_plot"):
    """
    Build a histogram to determine what the maxes look & visualize match_count
    Might be used to determine a resonable threshold

    @param coverage: the average coverage for an single nt
    @param show: Show visualization with match maxes
    @param save_fn: Save to disk with this file name or else it will be the default

    @return: the histogram array
    """
    #import matplotlib
    #matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    maxes = self.match_count.max(1) # get maxes along 1st dim

    h = plt.hist(maxes, bins=self.match_count.shape[0]) # figure out where the majority

    plt.ylabel("Frequency")
    plt.xlabel("Count per index")
    plt.title("Frequency count histogram")

    if show: plt.show()
    if save: plt.savefig(save_fn, dpi=160, frameon=False)

    return h[0]