def nova_plot():

	erg2mev=624151.

	fig=plot.figure()
	yrange = [1e-6,2e-4]
	xrange = [1e-1,1e5]
	plot.fill_between([0.2,10e3],[yrange[1],yrange[1]],[yrange[0],yrange[0]],facecolor='yellow',interpolate=True,color='yellow',alpha=0.5)
	plot.annotate('AMEGO',xy=(3,9e-5),xycoords='data',fontsize=26,color='black')

	lat=ascii.read("data/NMon2012.LAT.dat",names=['energy','en_low','en_high','flux','flux_err','tmp'])
	plot.scatter(lat['energy'],lat['flux']*erg2mev,color='red')
	plot.errorbar(lat['energy'],lat['flux']*erg2mev,xerr=[lat['en_low'],lat['en_high']],yerr=lat['flux_err']*erg2mev,ecolor='red',capsize=0,fmt='none')
	latul=ascii.read("data/NMon2012.LAT.limits.dat",names=['energy','en_low','en_high','flux','tmp1','tmp2','tmp3','tmp4'])
	plot.errorbar(latul['energy'],latul['flux']*erg2mev,xerr=[latul['en_low'],latul['en_high']],yerr=0.5*latul['flux']*erg2mev,uplims=True,ecolor='red',capsize=0,fmt='none')
	plot.scatter(latul['energy'],latul['flux']*erg2mev,color='red')

	leptonic=ascii.read("data/sp-NMon12-IC-best-fit-1MeV-30GeV.txt",names=['energy','flux'],data_start=1)
	hadronic=ascii.read("data/sp-NMon12-pi0-and-secondaries.txt",names=['energy','flux1','flux2'],data_start=1)	

	plot.plot(leptonic['energy'],leptonic['flux']*erg2mev,'r--',color='black',lw=2,label='Leptonic')
	plot.plot(hadronic['energy'],hadronic['flux2']*erg2mev,color='black',lw=2,label='Hadronic+Secondary Leptons')

	plot.legend(loc='upper right',fontsize='small',frameon=False,framealpha=0.5)
	plot.xscale('log')
	plot.yscale('log')
	plot.ylim(yrange)
	plot.xlim(xrange)
	plot.xlabel(r'Energy (MeV)')
	plot.ylabel(r'Energy$^2 \times $ Flux (Energy) (erg cm$^{-2}$ s$^{-1}$)')
	plot.title('Nova V339 Del 2013')
	plot.savefig('Nova_SED.png', bbox_inches='tight')
	plot.savefig('Nova_SED.eps', bbox_inches='tight')
	plot.show()
	plot.close()

def plot_prob_effector(sens, fpr, xmax=1, baserate=0.1):
    """Plots a line graph of P(effector|positive test) against
    the baserate of effectors in the input set to the classifier.
        
    The baserate argument draws an annotation arrow
    indicating P(pos|+ve) at that baserate
    """
    assert 0.1 <= xmax <= 1, "Max x axis value must be in range [0,1]"
    assert 0.01 <= baserate <= 1, "Baserate annotation must be in range [0,1]"
    baserates = pylab.arange(0, 1.05, xmax * 0.005)  
    probs = [p_correct_given_pos(sens, fpr, b) for b in baserates]
    pylab.plot(baserates, probs, 'r')
    pylab.title("P(eff|pos) vs baserate; sens: %.2f, fpr: %.2f" % (sens, fpr))
    pylab.ylabel("P(effector|positive)")
    pylab.xlabel("effector baserate")
    pylab.xlim(0, xmax)
    pylab.ylim(0, 1)
    # Add annotation arrow
    xpos, ypos = (baserate, p_correct_given_pos(sens, fpr, baserate))
    if baserate < xmax:
        if xpos > 0.7 * xmax:
            xtextpos = 0.05 * xmax
        else:
            xtextpos = xpos + (xmax-xpos)/5.
        if ypos > 0.5:
            ytextpos = ypos - 0.05
        else:
            ytextpos = ypos + 0.05
        pylab.annotate('baserate: %.2f, P(pos|+ve): %.3f' % (xpos, ypos), 
                       xy=(xpos, ypos), 
                       xytext=(xtextpos, ytextpos),
                       arrowprops=dict(facecolor='black', shrink=0.05))
    else:
        pylab.text(0.05 * xmax, 0.95, 'baserate: %.2f, P(pos|+ve): %.3f' %
                   (xpos, ypos))

def plot_trajectory(mu_vector):
    data0 = mu_vector[:, 0]
    data1 = mu_vector[:, 1]
    labels = ["{0}".format(i) for i in xrange(len(mu_vector))]
    plt.scatter(data0[:, 0], data0[:, 1], color="red")
    plt.scatter(data1[:, 0], data1[:, 1], color="blue")
    for i in xrange(len(mu_vector)):
        plt.annotate(
            labels[i],
            (data0[i, 0], data0[i, 1]),
            fontsize=5,
            xytext=(-10, 20),
            textcoords="offset points",
            ha="right",
            va="bottom",
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=0"),
        )
        plt.annotate(
            labels[i],
            (data1[i, 0], data1[i, 1]),
            fontsize=5,
            xytext=(-10, 20),
            textcoords="offset points",
            ha="right",
            va="bottom",
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=0"),
        )
    plt.savefig("Mean_Trajectory.png")
    plt.show()

    def test2DpyEI(self):
        
        f = lambda x: sum(sin(x))
        bounds = [[0., 5.], [0., 5.]]
        X = lhcSample(bounds, 5, seed=24)
        Y = [f(x) for x in X]

        kernel = GaussianKernel_ard(array([1.0, 1.0]))
        GP = GaussianProcess(kernel, X, Y)

        maxei = maximizeEI(GP, bounds)
        
        if False:
            figure(1)
            c0 = [(i/50.)*(bounds[0][1]-bounds[0][0])+bounds[0][0] for i in xrange(51)]
            c1 = [(i/50.)*(bounds[1][1]-bounds[1][0])+bounds[1][0] for i in xrange(51)]
            z = array([[GP.ei(array([i, j])) for i in c0] for j in c1])

            ax = plt.subplot(111)
            cs = ax.contour(c0, c1, z, 10, alpha=0.5, cmap=cm.Blues_r)
            plot([x[0] for x in X], [x[1] for x in X], 'ro')
            for i in xrange(len(X)):
                annotate('%2f'%Y[i], X[i])
            plot(maxei[1][0], maxei[1][1], 'ko')
            show()

def fancy_dendrogram(*args, **kwargs):
    '''
    Source: https://joernhees.de/blog/2015/08/26/scipy-hierarchical-clustering-and-dendrogram-tutorial/
    '''
    from scipy.cluster import hierarchy
    import matplotlib.pylab as plt
    
    max_d = kwargs.pop('max_d', None)
    if max_d and 'color_threshold' not in kwargs:
        kwargs['color_threshold'] = max_d
    annotate_above = kwargs.pop('annotate_above', 0)

    ddata = hierarchy.dendrogram(*args, **kwargs)

    if not kwargs.get('no_plot', False):
        plt.title('Hierarchical Clustering Dendrogram (truncated)')
        plt.xlabel('sample index or (cluster size)')
        plt.ylabel('distance')
        for i, d, c in zip(ddata['icoord'], ddata['dcoord'], ddata['color_list']):
            x = 0.5 * sum(i[1:3])
            y = d[1]
            if y > annotate_above:
                plt.plot(x, y, 'o', c=c)
                plt.annotate("%.3g" % y, (x, y), xytext=(0, -5),
                             textcoords='offset points',
                             va='top', ha='center')
        if max_d:
            plt.axhline(y=max_d, c='k')
    return ddata

def disk_plot(e, D, DrTh, color='g'):
    # plot disk to illustrate the weight strength
    if np.all(e==0):
        return
    # rescale to 0-1
    re = rescale(e) * D
    y,x = np.nonzero(re)
    r = re[(y,x)]
    # sort the disk from small to large
    locs = np.argsort( r )
    y = y[ locs ]
    x = x[ locs ]
    r = r[ locs ]
    # eleminate the small disks
    y = y[ r > DrTh ]
    x = x[ r > DrTh ]
    r = r[ r > DrTh ]
    plt.scatter(x,y,r, c=color, alpha=0.6, linewidths=0)

    # print the maximum value
    dev = max(x.max(), y.max()) * 0.07
    plt.annotate("%d" % e.max(), xy=(x[-1],y[-1]),
                xytext=(x[-1]+dev, y[-1]+dev),
                color = 'white',
                arrowprops=dict(color='white',
                arrowstyle="->"))

def measure_psf(vignet, pixscale=1., show=False, mask_value=None):
    y, x = np.mgrid[-vignet.shape[0]/2:vignet.shape[0]/2, -vignet.shape[1]/2:vignet.shape[1]/2]*pixscale
    if mask_value :
        vignet = ma.masked_values(vignet, mask_value).filled(0)
    # Fit the data using astropy.modeling
    p_init=models.Gaussian2D(amplitude=vignet.max(), x_mean=0., y_mean=0.,
        x_stddev=2*pixscale, y_stddev=2*pixscale, theta=0, cov_matrix=None)
    fit_p = fitting.LevMarLSQFitter()

    p = fit_p(p_init, x, y, vignet)
    barycenter=measure_barycenter(vignet, pixscale=pixscale)
    
    # Plot the data with the best-fit model
    P.figure(figsize=(8, 2.5))
    P.subplot(1, 3, 1)
    P.imshow(vignet, origin='lower', interpolation='nearest', vmin=vignet.min(), vmax=vignet.max())
    P.title("Data")
    P.subplot(1, 3, 2)
    P.imshow(p(x, y), origin='lower', interpolation='nearest', vmin=vignet.min(), vmax=vignet.max())
    P.scatter(vignet.shape[0]/2, vignet.shape[1]/2,marker="+")
    P.annotate("({:.3f},{:.3f})".format(*barycenter), (vignet.shape[0]/3, vignet.shape[1]/3))
    P.title("Model - psf = {:.2f}".format(2.3548*np.mean([p.x_stddev.value, p.y_stddev.value])))
    P.subplot(1, 3, 3)
    P.imshow(vignet - p(x, y), origin='lower', interpolation='nearest', vmin=-vignet.max()/10,vmax=vignet.max()/10)
    P.title("Residual")
    P.tight_layout()
    if show :
        P.show()
    
    return p

    def write_plots(self, take_every=1, show=False):
        pre = self.get_path_prefix()
        x, y, slope, yint = self.get_flux_xy(take_every)

        defaults = {'markersize':1, 'linewidth':0.1}
        fig = pl.figure()

        pl.plot(x, y, 'o',
                color='black', **defaults)

        fit_xs = np.array(pl.xlim())
        pl.plot(fit_xs, slope * fit_xs + yint,
                color='black', alpha=0.3, linewidth=1)

        pl.ylabel('Cumulative grains through aperture')
        pl.xlabel('Time')
        pl.title('Accumulation of grains vs time')

        pl.annotate('$Slope={:.1f}$'.format(slope),
            xy=(0.05, 0.95), xycoords='axes fraction', ha='left', va='top'
            , fontsize=14)
##        # Sim info:
##        pl.ylim(ymin=-2)  # make room for it
##        pl.annotate(info_for_naming, xy=(0.01, 0.03), xycoords='axes fraction', fontsize=12)
        pl.savefig('{}flux.png'.format(pre))
        pl.savefig('{}flux.svg'.format(pre))  # can take a while if complex plot
        if show:
            pl.show()
        else:
            pl.close()

def Plot2DQuadGeometry(xyz,IEN,nodeId,elemId):
    for i in np.arange(IEN.shape[0]):
        x = np.zeros(4)
        y = np.zeros(4)
        index = IEN[i,0]
        x[0] = xyz[index-1,0]
        y[0] = xyz[index-1,1]
        index = IEN[i,1]
        x[1] = xyz[index-1,0]
        y[1] = xyz[index-1,1] 
        index = IEN[i,2]
        x[2] = xyz[index-1,0]
        y[2] = xyz[index-1,1]
        x_mean = np.mean(x[0:3])
        y_mean = np.mean(y[0:3])
        x[3] = x[0]
        y[3] = y[0]
        plt.plot(x,y,color='black')
        if elemId == True:
            plt.annotate(str(i+1),(x_mean,y_mean))
    for i in np.arange(xyz.shape[0]):
        if nodeId == True:
            plt.plot(xyz[:,0],xyz[:,1],'go')
            plt.annotate(str(i+1),(xyz[i,0],xyz[i,1]))
            for i in np.arange(xyz.shape[0]):
                x = xyz[i,0]
                y = xyz[i,1]        

def plot(embeddings,labels):
    assert embeddings.shape[0] >= len(labels) , 'More labels than embeddings'
    pylab.figure(figsize=(15,15)) #in inches
    for i, label in enumerate(labels):
        x,y = embeddings[i,:]
        pylab.scatter(x,y)
        pylab.annotate(label,xy=(x,y),xytext=(5,2),textcoords='offset points',ha='right',va='bottom')
    pylab.show()

def plot(embeddings, labels, out):
    assert embeddings.shape[0] >= len(labels), 'More labels than embeddings'
    pylab.figure(figsize=(15,15))
    for i, label in enumerate(labels):
        x, y = embeddings[i,:]
        pylab.scatter(x, y)
        pylab.annotate(label, xy=(x, y), xytext=(5, 2), textcoords='offset points',
                       ha='right', va='bottom')
    pylab.savefig(out)
    print('Saved plot to {}'.format(out))

 def plot(self, basemap, annotate=True, lsize='small', kwargs_an=None, **kwargs_in):
     kwargs = dict(marker='o')
     kwargs.update(kwargs_in)
     if kwargs_an is None:
         kwargs_an = {}
     for key, val in self.items():
         x, y = basemap(val.longitude, val.latitude)
         basemap.plot((x,), (y,), **kwargs)
         if annotate:
             import matplotlib.pylab as plt
             plt.annotate(key, (x, y), xytext=(3, 3),
                          textcoords='offset points', size=lsize, **kwargs_an)

    def plot_frequency_altitude(self, f=2.0, ax=None, median_filter=False,
        vmin=None, vmax=None, altitude_range=(-99.9, 399.9), colorbar=False, return_image=False, annotate=True):

        if vmin is None:
            vmin = self.vmin

        if vmax is None:
            vmax = self.vmax

        if ax is None:
            ax = plt.gca()

        plt.sca(ax)
        plt.cla()
        freq_extent = (self.extent[0], self.extent[1],
            altitude_range[1], altitude_range[0])

        i = self.ionogram_list[0]
        inx = 1.0E6* (i.frequencies.shape[0] * f) / (i.frequencies[-1] - i.frequencies[0])
        img = self.tser_arr_all[:,int(inx),:]

        new_altitudes = np.arange(altitude_range[0], altitude_range[1], 14.)
        new_img = np.zeros((new_altitudes.shape[0], img.shape[1])) + np.nan

        for i in self.ionogram_list:
            e = int( round((i.time - self.extent[0]) / ais_code.ais_spacing_seconds ))

            pos = mex.iau_r_lat_lon_position(float(i.time))
            altitudes = pos[0] - ais_code.speed_of_light_kms * ais_code.ais_delays * 0.5 - mex.mars_mean_radius_km
            s = np.argsort(altitudes)
            new_img[:, e] = np.interp(new_altitudes, altitudes[s], img[s,e], left=np.nan, right=np.nan)

        plt.imshow(new_img, vmin=vmin, vmax=vmax,
            interpolation='Nearest', extent=freq_extent, origin='upper', aspect='auto')

        plt.xlim(freq_extent[0], freq_extent[1])
        plt.ylim(*altitude_range)

        ax.set_xlim(self.extent[0], self.extent[1])
        ax.xaxis.set_major_locator(celsius.SpiceetLocator())

        celsius.ylabel(r'Alt./km')
        if annotate:
            plt.annotate('f = %.1f MHz' % f, (0.02, 0.9),
                    xycoords='axes fraction', color='cyan', verticalalignment='top', fontsize='small')

        if colorbar:
            old_ax = plt.gca()
            plt.colorbar(cax = celsius.make_colorbar_cax(), ticks=self.cbar_ticks).set_label(r"$Log_{10} V^2 m^{-2} Hz^{-1}$")
            plt.sca(old_ax)

        if return_image:
            return new_img, freq_extent, new_altitudes

def plot(words):
    %matplotlib inline
    embeddings = [model[w] for w in words]

    pca = PCA(n_components=2)  
    two_d_embeddings = pca.fit_transform(embeddings)

    pylab.figure(figsize=(5,5))  # in inches
    for i, label in enumerate(words):
        x, y = two_d_embeddings[i,:]
        pylab.scatter(x, y)
        pylab.annotate(label, xy=(x, y), xytext=(5, 2), textcoords='offset points', ha='right', va='bottom')
    pylab.show()

def plot_embedding(embeddings, labels):
    """Applies non-linear dimensionalty reduction using tSNE and plots
    the words."""
    tsne = TSNE(perplexity=30, n_components=2, init='pca', n_iter=5000)
    two_d_embeddings = tsne.fit_transform(embeddings)

    pylab.figure(figsize=(15,15))  # in inches
    for i, label in enumerate(labels):
        x, y = two_d_embeddings[i,:]
        pylab.scatter(x, y)
        pylab.annotate(label, xy=(x, y), xytext=(5, 2),
            textcoords='offset points', ha='right', va='bottom')
    pylab.show()

def plotSamples(samples, marker):
    '''
    绘制样本
    '''
    xVals = []
    yVals = []
    for s in samples:
        x = s.getFeatures()[0]
        y = s.getFeatures()[1]
        plt.annotate(s.getName(), xy=(x,y), xytext=(x+0.13, y-0.07), fontsize='x-large')
        xVals.append(x)
        yVals.append(y)
    plt.plot(xVals, yVals, marker)

def plot(words):
    embeddings = [model[w] for w in words]

    pca = PCA(n_components = 2)
    two_d_embeddings = pca.fit_transform(embeddings)

    pylab.figure(figsize=(5,5))
    for i, label in enumerate(words):
        x, y = two_d_embeddings[i,:]
        pylab.scatter(x,y)
        pylab.annotate(label, xy=(x,y), xytext=(28, 2), 
            textcoords = 'offset points', ha='right', va='bottom')
    pylab.savefig('testplot19.png')

def check_psp_fitting():
    """Plots the results of the current fitting with the save fits and denotes 
    when there is a change. 
    """
    plotting=True # specifies whether to make plots of fitting results
        
    test_data_files=[os.path.join(test_data_dir,f) for f in os.listdir(test_data_dir)] #list of test files
    for file in sorted(test_data_files):
#    for file in ['test_psp_fit/1492546902.92_2_6stacked.json']: order of parameters affects this fit
        print 'file', file
        test_dict=json.load(open(file)) # load test data
        avg_trace=neuroanalysis.data.Trace(data=np.array(test_dict['input']['data']), dt=test_dict['input']['dt']) # create Trace object
        psp_fits = fit_psp(avg_trace, 
                           sign=test_dict['input']['amp_sign'], 
                           stacked=test_dict['input']['stacked'] 
                            )                        
        
        change_flag=False
        if test_dict['out']['best_values']!=psp_fits.best_values:     
            print '  the best values dont match'
            print '\tsaved', test_dict['out']['best_values']
            print '\tobtained', psp_fits.best_values
            change_flag=True
            
        if test_dict['out']['best_fit']!=psp_fits.best_fit.tolist():
            print '  the best fit traces dont match'
            print '\tsaved', test_dict['out']['best_fit']
            print '\tobtained', psp_fits.best_fit.tolist()
            change_flag=True
        
        if test_dict['out']['nrmse']!=float(psp_fits.nrmse()):
            print '  the nrmse doesnt match'
            print '\tsaved', test_dict['out']['nrmse']
            print '\tobtained', float(psp_fits.nrmse())
            change_flag=True
            
        if plotting:
            import matplotlib.pylab as mplt
            fig=mplt.figure(figsize=(20,8))
            ax=fig.add_subplot(1,1,1)
            ax2=ax.twinx()
            ax.plot(avg_trace.time_values, psp_fits.data*1.e3, 'b', label='data')
            ax.plot(avg_trace.time_values, psp_fits.best_fit*1.e3, 'g', lw=5, label='current best fit')
            ax2.plot(avg_trace.time_values, test_dict['input']['weight'], 'r', label='weighting')
            if change_flag is True:
                ax.plot(avg_trace.time_values, np.array(test_dict['out']['best_fit'])*1.e3, 'k--', lw=5, label='original best fit')
                mplt.annotate('CHANGE', xy=(.5, .5), xycoords='figure fraction', fontsize=40)
            ax.legend()
            mplt.title(file + ', nrmse =' + str(psp_fits.nrmse()))
            mplt.show()

def pkmean(kmean, X, label_centers=True, save_loc=False):
    '''
    INPUT: fitted kmean model,
           1-dimensional array of x,y tuples,
           True/False to include centroid labels,
           save location
    OUTPUT: None

    Create and show a scatterplot color points by cluster
    and plot the cluster centers in black.

    You can choose to save the plot if
    file location is specified for save_loc
    '''
    centers = kmean.cluster_centers_
    colors = plt.cm.Spectral(np.linspace(0, 1, len(centers)))
    labels = np.arange(len(centers))
    # plot and color points by cluster
    for label, col in zip(labels, colors):
        label_indices = np.where(kmean.labels_ == label)
        x = X[:, 0][label_indices]
        y = X[:, 1][label_indices]

        plt.plot(x, y, '.',
                 markerfacecolor=col,
                 alpha=0.1)
        plt.axis('off')

    # add labels to centers
    if label_centers:
        center_x = centers[:, 0]
        center_y = centers[:, 1]
        # plot cluster centers in black
        plt.plot(center_x, center_y, '.',
                 markerfacecolor='k')

        for label, x, y in zip(labels, center_x, center_y):
            plt.annotate(
                label,
                xy=(x, y), xytext = (-20, 20),
                textcoords = 'offset points', ha = 'right', va = 'bottom',
                bbox = dict(boxstyle='round,pad=0.5', fc='yellow', alpha=0.5),
                arrowprops = dict(arrowstyle='->',
                                  connectionstyle='arc3,rad=0'))

    if save_loc:
        plt.savefig(save_loc)
    else:
        plt.show()

def plot_word_embedding(plt, table, labels=None, title='', num=1):
    import numpy as np
    plt.figure(num)
    vectors = np.matrix(table).tolist()
    words = list(table.index)

    import matplotlib
    if(type(labels) == type(None)):
        None
        colors = None
    else:
        label_set = list(set(list(labels.values.transpose().tolist())[0]))

        def get_spaced_colors(n):
            max_value = 16581375  # 255**3
            interval = int(max_value / n)
            colors = [hex(I)[2:].zfill(6)
                      for I in range(0, max_value, interval)]

            return [(int(i[:2], 16), int(i[2:4], 16), int(i[4:], 16)) for i in colors]
        colors = get_spaced_colors(len(label_set))

    for i in range(len(words)):
        point = vectors[i]
        word = words[i]
        # plot points
        plt.scatter(point[0], point[1])
        # plot word annotations
        if(type(labels) == type(None)):

            plt.annotate(
                word,
                xy=(point[0], point[1]),
                size="x-small"
            )
        else:
            label_index = label_set.index(
                list(labels.values.transpose().tolist())[0][i])
            plt.annotate(
                word,
                xy=(point[0], point[1]),
                color='#' +
                "".join(list(map(lambda x: format(x, '#04x')
                                 [2:], colors[label_index]))).upper(),
                size="x-small"
            )

    plt.tight_layout()
    plt.title(title)