def plot_call_rate(c):
    # Histogram
    P.clf()
    P.figure(1)
    P.hist(c[:,1], normed=True)
    P.xlabel('Call Rate')
    P.ylabel('Portion of Variants')
    P.savefig(os.environ['OBER'] + '/doc/imputation/cgi/call_rate.png')

####################################################################################
#if __name__ == '__main__':
#    # Input parameters
#    file_name = sys.argv[1]  # Name of data file with MAF, call rates
#
#    # Load data
#    c = np.loadtxt(file_name, dtype=np.float16)
#
#    # Breakdown by call rate (proportional to the #samples, 1415)
#    plot_call_rate(c)
#    h = np.histogram(c[:,1])
#    a = np.flipud(np.cumsum(np.flipud(h[0])))/float(c.shape[0])
#    print np.concatenate((h[1][:-1][newaxis].transpose(), a[newaxis].transpose()), axis=1)

    # Breakdown by minor allele frequency
    maf_n = 20
    maf_bins = np.linspace(0, 0.5, maf_n + 1)
    maf_bin = np.digitize(c[:,0], maf_bins)
    d = c.astype(float64)
    mean_call_rate = np.array([(1.*np.mean(d[maf_bin == i,1])) for i in xrange(len(maf_bins))])
    P.bar(maf_bins - h, mean_call_rate, width=h)

    P.figure(2)
    h = (maf_bins[-1] - maf_bins[0]) / maf_n
    P.bar(maf_bins - h, mean_call_rate, width=h)
    P.savefig(os.environ['OBER'] + '/doc/imputation/cgi/call_rate_maf.png')

def study_multiband_planck(quick=True):
    savename = datadir+'cl_multiband.pkl'
    bands = [100, 143, 217, 'mb']
    if quick: cl = pickle.load(open(savename,'r'))
    else:
        cl = {}
        mask = load_planck_mask()
        mask_factor = np.mean(mask**2.)
        for band in bands:
            this_map = load_planck_data(band)
            this_cl = hp.anafast(this_map*mask, lmax=lmax)/mask_factor
            cl[band] = this_cl
        pickle.dump(cl, open(savename,'w'))


    cl_theory = {}
    pl.clf()
    
    for band in bands:
        l_theory, cl_theory[band] = get_cl_theory(band)
        this_cl = cl[band]
        pl.plot(this_cl/cl_theory[band])
        
    pl.legend(bands)
    pl.plot([0,4000],[1,1],'k--')
    pl.ylim(.7,1.3)
    pl.ylabel('data/theory')

def plot_many_corr_delta_vel():
    pl.clf()
    leg = []
    for kmax in [0.05, 0.1, 0.2, 0.5, 1., 2., 5.]:
        plot_corr_delta_vel(kmin=1e-3, kmax=kmax, doclf=False)
        leg.append('kmax=%0.2f'%kmax)
    pl.legend(leg)

def plot_mock(mock):
    plt.clf()
    plt.plot(mock['dates'], mock['y'], marker='+', color='blue',
             label='data', markersize=9)
    plt.plot(mock['dates'], mock['y_without_seasonal'],
             color='green', alpha=0.6, linewidth=1,
             label='model without seasonal')

def study_redmapper_lrg_3d(hemi='north'):
    # create 3d grid object
    grid = grid3d(hemi=hemi)
    
    # load SDSS data
    sdss = load_sdss_data_both_catalogs(hemi)
    
    # load redmapper catalog
    rm = load_redmapper(hemi=hemi)
    
    # get XYZ positions (Mpc) of both datasets
    x_sdss, y_sdss, z_sdss = grid.xyz_from_radecz(sdss['ra'], sdss['dec'], sdss['z'], applyzcut=False)
    x_rm, y_rm, z_rm = grid.xyz_from_radecz(rm['ra'], rm['dec'], rm['z_spec'], applyzcut=False)
    pos_sdss = np.vstack([x_sdss, y_sdss, z_sdss]).T
    pos_rm = np.vstack([x_rm, y_rm, z_rm]).T

    # build a couple of KDTree's, one for SDSS, one for RM.
    from sklearn.neighbors import KDTree
    tree_sdss = KDTree(pos_sdss, leaf_size=30)
    tree_rm = KDTree(pos_rm, leaf_size=30)

    lrg_counts = tree_sdss.query_radius(pos_rm, 100., count_only=True)
    pl.clf()
    pl.hist(lrg_counts, bins=50)
    
    
    ipdb.set_trace()

    def plot(self, bit_stream):
        if self.previous_bit_stream != bit_stream.to_list():
            self.previous_bit_stream = bit_stream

            x = []
            y = []
            bit = None

            for bit_time in bit_stream.to_list():
                if bit is None:
                    x.append(bit_time)
                    y.append(0)
                    bit = 0
                elif bit == 0:
                    x.extend([bit_time, bit_time])
                    y.extend([0, 1])
                    bit = 1
                elif bit == 1:
                    x.extend([bit_time, bit_time])
                    y.extend([1, 0])
                    bit = 0

            plt.clf()
            plt.plot(x, y)
            plt.xlim([0, 10000])
            plt.ylim([-0.1, 1.1])
            plt.show()
            plt.pause(0.005)

def PlotTurbulenceIllustr(a):

    """
    Can generate the grid with
    g=kolmogorovutils.GenerateKolmogorov3D( 1025, 129, 129)
    a=kolmogorovutils.GridToNumarray(g)

    """

    for x in [1,10,100]:

        suba= numarray.sum(a[:,:,0:x], axis=2)

        suba.transpose()
        pylab.clf()
        pylab.matshow(suba)
        pylab.savefig("temp/turb3d-sum%03i.eps" % x)

    for x in [1,10,100]:

        for j in [0,1,2]:

            suba= numarray.sum(a[:,:200,j*x:(j+1)*x], axis=2)

            suba.transpose()
            pylab.clf()
            pylab.matshow(suba)
            pylab.savefig("temp/turb3d-sum%03i-s%i.eps" % (x,j))        

def plot_average(filenames, save_plot=True, show_plot=False, dpi=100):

    ''' Plot Signal average from a list of averaged files. '''

    fname = get_files_from_list(filenames)

    # plot averages
    pl.ioff()  # switch off (interactive) plot visualisation
    factor = 1e15
    for fnavg in fname:
        name = fnavg[0:len(fnavg) - 4]
        basename = os.path.splitext(os.path.basename(name))[0]
        print fnavg
        # mne.read_evokeds provides a list or a single evoked based on condition.
        # here we assume only one evoked is returned (requires further handling)
        avg = mne.read_evokeds(fnavg)[0]
        ymin, ymax = avg.data.min(), avg.data.max()
        ymin *= factor * 1.1
        ymax *= factor * 1.1
        fig = pl.figure(basename, figsize=(10, 8), dpi=100)
        pl.clf()
        pl.ylim([ymin, ymax])
        pl.xlim([avg.times.min(), avg.times.max()])
        pl.plot(avg.times, avg.data.T * factor, color='black')
        pl.title(basename)

        # save figure
        fnfig = os.path.splitext(fnavg)[0] + '.png'
        pl.savefig(fnfig, dpi=dpi)

    pl.ion()  # switch on (interactive) plot visualisation

    def plotFittingResults(self):
        """
        Plot results of Rmax optimization procedure and best fit of the experimental data
        """
        _listFitQ = [tmp.getValue() for tmp in self.getDataOutput().getScatteringFitQ()]
        _listFitValues = [tmp.getValue() for tmp in self.getDataOutput().getScatteringFitValues()]
        _listExpQ = [tmp.getValue() for tmp in self.getDataInput().getExperimentalDataQ()]
        _listExpValues = [tmp.getValue() for tmp in self.getDataInput().getExperimentalDataValues()]

        #_listExpStdDev = None
        #if self.getDataInput().getExperimentalDataStdDev():
        #    _listExpStdDev = [tmp.getValue() for tmp in self.getDataInput().getExperimentalDataStdDev()]
        #if _listExpStdDev:
        #    pylab.errorbar(_listExpQ, _listExpValues, yerr=_listExpStdDev, linestyle='None', marker='o', markersize=1,  label="Experimental Data")
        #    pylab.gca().set_yscale("log", nonposy='clip')
        #else:         
        #    pylab.semilogy(_listExpQ, _listExpValues, linestyle='None', marker='o', markersize=5,  label="Experimental Data")

        pylab.semilogy(_listExpQ, _listExpValues, linestyle='None', marker='o', markersize=5, label="Experimental Data")
        pylab.semilogy(_listFitQ, _listFitValues, label="Fitting curve")
        pylab.xlabel('q')
        pylab.ylabel('I(q)')
        pylab.suptitle("RMax : %3.2f. Fit quality : %1.3f" % (self.getDataInput().getRMax().getValue(), self.getDataOutput().getFitQuality().getValue()))
        pylab.legend()
        pylab.savefig(os.path.join(self.getWorkingDirectory(), "gnomFittingResults.png"))
        pylab.clf()

 def plot(self,key='Re'):
     """
     Create a plot of a variable over the ORACLES study area. 
     
     Parameters
     ----------
     key : string
     See names for available datasets to plot.
     
     clf : boolean
     If True, clear off pre-existing figure. If False, plot over pre-existing figure.
     
     Modification history
     --------------------
     Written: Michael Diamond, 08/16/2016, Seattle, WA
     Modified: Michael Diamond, 08/21/2016, Seattle, WA
        -Added ORACLES routine flight plan, Walvis Bay (orange), and Ascension Island
     Modified: Michael Diamond, 09/02/2016, Swakopmund, Namibia
         -Updated flihgt track
     """
     plt.clf()
     size = 16
     font = 'Arial'
     m = Basemap(llcrnrlon=self.lon.min(),llcrnrlat=self.lat.min(),urcrnrlon=self.lon.max(),\
     urcrnrlat=self.lat.max(),projection='merc',resolution='i')
     m.drawparallels(np.arange(-180,180,5),labels=[1,0,0,0],fontsize=size,fontname=font)
     m.drawmeridians(np.arange(0,360,5),labels=[1,1,0,1],fontsize=size,fontname=font)
     m.drawmapboundary(linewidth=1.5)        
     m.drawcoastlines()
     m.drawcountries()
     if key == 'Pbot' or key == 'Ptop' or key == 'Nd' or key == 'DZ': 
         m.drawmapboundary(fill_color='steelblue')
         m.fillcontinents(color='floralwhite',lake_color='steelblue',zorder=0)
     else: m.fillcontinents('k',zorder=0)
     if key == 'Nd':
         m.pcolormesh(self.lon,self.lat,self.ds['%s' % key],cmap=self.colors['%s' % key],\
         latlon=True,norm = LogNorm(vmin=self.v['%s' % key][0],vmax=self.v['%s' % key][1]))
     elif key == 'Zbf' or key == 'Ztf':
         levels = [0,250,500,750,1000,1250,1500,1750,2000,2500,3000,3500,4000,5000,6000,7000,8000,9000,10000]
         m.contourf(self.lon,self.lat,self.ds['%s' % key],levels=levels,\
         cmap=self.colors['%s' % key],latlon=True,extend='max')
     elif key == 'DZ':
         levels = [0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000]
         m.contourf(self.lon,self.lat,self.ds['%s' % key],levels=levels,\
         cmap=self.colors['%s' % key],latlon=True,extend='max')
     else:
         m.pcolormesh(self.lon,self.lat,self.ds['%s' % key],cmap=self.colors['%s' % key],\
         latlon=True,vmin=self.v['%s' % key][0],vmax=self.v['%s' % key][1])
     cbar = m.colorbar()
     cbar.ax.tick_params(labelsize=size-2) 
     cbar.set_label('[%s]' % self.units['%s' % key],fontsize=size,fontname=font)
     if key == 'Pbot' or key == 'Ptop': cbar.ax.invert_yaxis() 
     m.scatter(14.5247,-22.9390,s=250,c='orange',marker='D',latlon=True)
     m.scatter(-14.3559,-7.9467,s=375,c='c',marker='*',latlon=True)
     m.scatter(-5.7089,-15.9650,s=375,c='chartreuse',marker='*',latlon=True)
     m.plot([14.5247,13,0],[-22.9390,-23,-10],c='w',linewidth=5,linestyle='dashed',latlon=True)
     m.plot([14.5247,13,0],[-22.9390,-23,-10],c='k',linewidth=3,linestyle='dashed',latlon=True)
     plt.title('%s from MSG SEVIRI on %s/%s/%s at %s UTC' % \
     (self.names['%s' % key],self.month,self.day,self.year,self.time),fontsize=size+4,fontname=font)
     plt.show()

def plotLSQ(C_ms,lsqSc,lsqMu,lsqSl,LSQ,viewDirectory,TextSize=16):
  C_MS = C_ms - 1.0
  #480x480
  myfile = os.path.join(viewDirectory,'lsq.png')
  title='Cumulative LSQ Penalties' 
  red = 'S(Q,E)'
  green = 'high E scatter'
  blue = 'slope of high E scatter'
  xlabel = r"$C_{ms}$ [unitless]"
  ylabel = "Cumulative Penalty [unitless]"
  pylab.clf()
  pylab.plot( C_MS, lsqSc, 'r-', linewidth=5 )
  pylab.plot( C_MS, lsqMu, 'g-', linewidth=5 )
  pylab.plot( C_MS, lsqSl, 'b-', linewidth=5 )
  pylab.grid( 1 )
  pylab.legend( (red, green, blue), loc="upper left" )
  pylab.xlabel( xlabel )
  pylab.ylabel( ylabel )
  pylab.title(title)
  pylab.savefig( myfile )
  
  myfile = os.path.join(viewDirectory,'lsq_f.png')
  pylab.clf()
  title='Final Cumulative LSQ Penalty'
  pylab.plot( C_MS,LSQ, 'b-', linewidth = 5)
  pylab.xlabel( xlabel )
  pylab.ylabel( ylabel )
  pylab.title( title )
  pylab.grid(1)
  pylab.savefig( myfile )
  return

def plotRocCurves(lesion, lesion_en):
	file_legend = []
	for techniqueMid in techniquesMid:
		for techniqueLow in techniquesLow:
			file_legend.append((directory + techniqueLow + "/" + techniqueMid + "/operating-points-" + lesion + "-scale.dat", "Low-level: " + techniqueLow + ". Mid-level: " + techniqueMid + "."))
			
			pylab.clf()
			pylab.figure(1)
			pylab.xlabel('1 - Specificity', fontsize=12)
			pylab.ylabel('Sensitivity', fontsize=12)
			pylab.title(lesion_en)
			pylab.grid(True, which='both')
			pylab.xticks([i/10.0 for i in range(1,11)])
			pylab.yticks([i/10.0 for i in range(0,11)])
			#pylab.tick_params(axis="both", labelsize=15)
			
			for file, legend in file_legend:
				points = open(file,"rb").readlines()
				x = [float(p.split()[0]) for p in points]
				y = [float(p.split()[1]) for p in points]
				x.append(0.0)
				y.append(0.0)
				
				auc = numpy.trapz(y, x) * -100

				pylab.grid()
				pylab.plot(x, y, '-', linewidth = 1.5, label = legend + u" (AUC = {0:0.1f}%)".format(auc))

	pylab.legend(loc = 4, borderaxespad=0.4, prop={'size':12})
	pylab.savefig(directory + "plots/" + lesion + ".pdf", format='pdf')

    def plot_df(self,show=False):
        from matplotlib import pylab as plt 

        if self.afp is None:
            print 'afp not initilized. call update afp'
            return -1

        linecords,td,df,rtn,minmaxy = self.afp

        formatter = PlotDateFormatter(df.index)
        #fig = plt.figure()
        #ax = plt.addsubplot()
        fig, ax = plt.subplots()
        ax.xaxis.set_major_formatter(formatter)
    
        ax.plot(np.arange(len(df)), df['p'])

        for cord in linecords:
            plt.plot(cord[0],cord[1],color='red')

        fig.autofmt_xdate()
        plt.xlim(-10,len(df.index) + 10)
        plt.ylim(df.p.min() - 10,df.p.max() + 10)
        plt.grid(ax)
        #if show:
        #    plt.show()
        
        #"{0}{1}.png".format("./data/",datetime.datetime.strftime(datetime.datetime.now(),'%Y%M%m%S'))
        if self.plot_file:
            save_path = self.plot_file.format(self.symbol)
            if os.path.exists(os.path.dirname(save_path)):
                plt.savefig(save_path)
                
        plt.clf()
        plt.close()

def check_HDF5(size=64):
    """
    Plot images with landmarks to check the processing
    """

    # Get hdf5 file
    hdf5_file = os.path.join(data_dir, "CelebA_%s_data.h5" % size)

    with h5py.File(hdf5_file, "r") as hf:
        data_color = hf["training_color_data"]
        data_lab = hf["training_lab_data"]
        data_black = hf["training_black_data"]
        for i in range(data_color.shape[0]):
            fig = plt.figure()
            gs = gridspec.GridSpec(3, 1)
            for k in range(3):
                ax = plt.subplot(gs[k])
                if k == 0:
                    img = data_color[i, :, :, :].transpose(1,2,0)
                    ax.imshow(img)
                elif k == 1:
                    img = data_lab[i, :, :, :].transpose(1,2,0)
                    img = color.lab2rgb(img)
                    ax.imshow(img)
                elif k == 2:
                    img = data_black[i, 0, :, :] / 255.
                    ax.imshow(img, cmap="gray")
            gs.tight_layout(fig)
            plt.show()
            plt.clf()
            plt.close()

def plotMassFunction(im, pm, outbase, mmin=9, mmax=13, mstep=0.05):
    """
    Make a comparison plot between the input mass function and the 
    predicted projected correlation function
    """
    plt.clf()

    nmbins = ( mmax - mmin ) / mstep
    mbins = np.logspace( mmin, mmax, nmbins )
    mcen = ( mbins[:-1] + mbins[1:] ) /2
    
    plt.xscale( 'log', nonposx = 'clip' )
    plt.yscale( 'log', nonposy = 'clip' )
    
    ic, e, p = plt.hist( im, mbins, label='Original Halos', alpha=0.5, normed = True)
    pc, e, p = plt.hist( pm, mbins, label='Added Halos', alpha=0.5, normed = True)
    
    plt.legend()
    plt.xlabel( r'$M_{vir}$' )
    plt.ylabel( r'$\frac{dN}{dM}$' )
    #plt.tight_layout()
    plt.savefig( outbase+'_mfcn.png' )
    
    mdtype = np.dtype( [ ('mcen', float), ('imcounts', float), ('pmcounts', float) ] )
    mf = np.ndarray( len(mcen), dtype = mdtype )
    mf[ 'mcen' ] = mcen
    mf[ 'imcounts' ] = ic
    mf[ 'pmcounts' ] = pc

    fitsio.write( outbase+'_mfcn.fit', mf )

 def psd(self, lc_data, n=1024, save=False):
     flux = lc_data[:,1]
     # check NaN value, and reset it to zero
     nani = np.where(np.isnan(flux) > 0)
     flux[nani] = 0.0
     
     thz = np.fft.fftfreq(n, d=self.exposuretime)
     fre = []
     for index in range(len(flux)/n):
         temp = np.fft.fft(flux[index*n:(index + 1)*n])
         fre.append(temp[1:n/2])
     fre_array = np.asarray(fre)
     avg_fre_array = np.average(fre_array, axis=0)
     power = (np.abs(avg_fre_array)**2)/thz[1:n/2]
     
     plt.clf()
     ps = plt.gca()    
     ps.plot(thz[1:n/2], power, 'k.-')
     ps.set_xscale('log')
     ps.set_yscale('log')
     ps.set_xlabel('Frequency (Hz)')
     ps.set_ylabel(r'PSD (power$^2$/frequency)')
     ps.set_title('PSD: ' + op.basename(self.filename) + r', n=%4d, $\Delta$t=%4.2f sec' % (n, self.exposuretime))
     plt.show()
     
     if save == 1:
         save_filename = 'f' + str(self.fibre) + '_' + op.basename(self.filename)[:-5] + '_psd.png'
         plt.savefig(save_filename, format='png')

def get_histogram_scale(distances_dict, nbins):
    """Draws histogram to outfile_name.
    """
    scale_dict = defaultdict(list)
    #draw histograms
    for d_dict in distances_dict.values():
        for i, (field, data) in enumerate(d_dict.items()):
            if len(data) < 1:
                continue
            histogram = hist(data,bins=nbins)
            
            fig  = gcf()
            axis = fig.gca()

            #get height scale: y/x
            ymin,ymax = axis.get_ylim()
        
            xmin,xmax = axis.get_xlim()
            scale_dict['ymin'].append(ymin)
            scale_dict['ymax'].append(ymax)
            scale_dict['xmin'].append(xmin)
            scale_dict['xmax'].append(xmax)

            clf()
    
    yscale = (min(scale_dict['ymin']),max(scale_dict['ymax']))
    xscale = (min(scale_dict['xmin']),max(scale_dict['xmax']))
    
    return xscale,yscale

def visualize_representations(dir, nlayers):
    pylab.clf()

    n_subplot_vertical = nlayers + 1
    n_subplot_horizontal = 4

    # input
    location = 1 + nlayers * n_subplot_horizontal + 1
    filename = dir + "representations/input.txt"
    plot_representation(n_subplot_vertical, n_subplot_horizontal, location, filename, "x")

    #
    for i in range(nlayers):
        # reconstruction of input
        location = 1 + (nlayers - i) * n_subplot_horizontal
        filename = dir + "representations/recons_from_hidden_l" + str(i) + ".txt"
        plot_representation(n_subplot_vertical, n_subplot_horizontal, location, filename, "rebuilt from h")

        # hidden layer
        location = 1 + (nlayers - i - 1) * n_subplot_horizontal + 1
        filename = dir + "representations/hidden_l" + str(i) + ".txt"
        plot_representation(n_subplot_vertical, n_subplot_horizontal, location, filename, "hidden")

        # reconstruction of layer through speech
        location = 1 + (nlayers - i - 1) * n_subplot_horizontal + 2
        filename = dir + "representations/recons_from_speech_l" + str(i) + ".txt"
        plot_representation(n_subplot_vertical, n_subplot_horizontal, location, filename, "rebuilt from s")

        # speech
        location = 1 + (nlayers - i - 1) * n_subplot_horizontal + 3
        filename = dir + "representations/speech_l" + str(i) + ".txt"
        plot_representation(n_subplot_vertical, n_subplot_horizontal, location, filename, "speech")

    pylab.savefig(dir + "representations.png")

def plot_feat_hist(data_name_list, filename=None):
    pylab.clf()
    num_rows = 1 + (len(data_name_list) - 1) / 2
    num_cols = 1 if len(data_name_list) == 1 else 2
    pylab.figure(figsize=(5 * num_cols, 4 * num_rows))

    for i in range(num_rows):
        for j in range(num_cols):
            pylab.subplot(num_rows, num_cols, 1 + i * num_cols + j)
            x, name = data_name_list[i * num_cols + j]
            pylab.title(name)
            pylab.xlabel('Value')
            pylab.ylabel('Density')
            # the histogram of the data
            max_val = np.max(x)
            if max_val <= 1.0:
                bins = 50
            elif max_val > 50:
                bins = 50
            else:
                bins = max_val
            n, bins, patches = pylab.hist(
                x, bins=bins, normed=1, facecolor='green', alpha=0.75)

            pylab.grid(True)

    if not filename:
        filename = "feat_hist_%s.png" % name

    pylab.savefig(os.path.join(CHART_DIR, filename), bbox_inches="tight")

def plotRocCurves(file_legend):
	pylab.clf()
	pylab.figure(1)
	pylab.xlabel('1 - Specificity', fontsize=12)
	pylab.ylabel('Sensitivity', fontsize=12)
	pylab.title("Need for Referral")
	pylab.grid(True, which='both')
	pylab.xticks([i/10.0 for i in range(1,11)])
	pylab.yticks([i/10.0 for i in range(0,11)])
	pylab.tick_params(axis="both", labelsize=15)

	for file, legend in file_legend:
		points = open(file,"rb").readlines()
		x = [float(p.split()[0]) for p in points]
		y = [float(p.split()[1]) for p in points]
		dev = [float(p.split()[2]) for p in points]
		x = [0.0] + x
		y = [0.0] + y
		dev = [0.0] + dev
	
		auc = np.trapz(y, x) * 100
		aucDev = np.trapz(dev, x) * 100

		pylab.grid()
		pylab.errorbar(x, y, yerr = dev, fmt='-')
		pylab.plot(x, y, '-', linewidth = 1.5, label = legend + u" (AUC = {0:0.1f}% \xb1 {1:0.1f}%)".format(auc,aucDev))

	pylab.legend(loc = 4, borderaxespad=0.4, prop={'size':12})
	pylab.savefig("referral/referral-curves.pdf", format='pdf')