def infer_latent_dim(X, verbose=0, maxr=-1):
    """
    r = infer_latent_dim(X, verbose=0)
    Infer the latent dimension of an aray assuming data+gaussian noise mixture
    
    Parameters
    ----------
    array X, data whose deimsnionhas to be inferred
    verbose=0, int, verbosity level
    maxr=-1, int, maximum dimension that can be achieved
             if maxr = -1, this is equal to rank(X)
    
    Returns
    -------
    r, int, the inferred dimension
    """
    if maxr ==-1:
        maxr = np.minimum(X.shape[0],X.shape[1])
        
    U,S,V = nl.svd(X,0)
    if verbose>1:
        print "Singular Values", S
    L = []
    for k in range(maxr):
        L.append(_linear_dim_criterion_(S,k,X.shape[1],X.shape[0])/X.shape[0])

    L = np.array(L)
    rank = np.argmax(L)
    if verbose:
        import matplotlib.pylab as mp
        mp.figure()
        mp.bar(np.arange(maxr),L-L.mean())

    return rank

def stackedHistogram(data, bins=10, range=None, log = False, normed = False,
                     filename = None, labels = None,
                     **formating):

    histograms = []
    d1 = data[0]
    curHist, bins = numpy.histogram(d1, bins = bins, range = range, normed = normed)
    histograms.append(curHist)

    for d in data[1:]:
        curHist, bins = numpy.histogram(d, bins = bins, normed = normed)
        histograms.append(curHist)

    width = bins[1] - bins[0]

    colors = 'b r k g c m y w'.split()
    nRepeats = len(data) / len(colors)
    for i in xrange(nRepeats):
        colors = colors.extend(colors)

    for histo, color in zip(histograms, colors):

        pylab.bar(bins[:-1], histo[:-1], width = width, log=log,
                  edgecolor = color, facecolor = None)


        
        
    doFormating(**formating)
    pylab.show()
    if filename is not None:
        pylab.savefig(filename)
        pylab.clf()

def plot_tuning_curves(direction_rates, title):
    """
    This function takes the x-values and the y-values  in units of spikes/s 
    (found in the two columns of direction_rates) and plots a histogram and 
    polar representation of the tuning curve. It adds the given title.
    """
    x = direction_rates[:,0]
    y = direction_rates[:,1]
    plt.figure()
    plt.subplot(2,2,1)
    plt.bar(x,y,width=45,align='center')
    plt.xlim(-22.5,337.5)
    plt.xticks(x)
    plt.xlabel('Direction of Motion (degrees)')
    plt.ylabel('Firing Rate (spikes/s)')
    plt.title(title)   
        
        
    
    plt.subplot(2,2,2,polar=True)
    r = np.append(y,y[0])
    theta = np.deg2rad(np.append(x, x[0]))
    plt.polar(theta,r,label='Firing Rate (spikes/s)')
    plt.legend(loc=8)
    plt.title(title)

    def eqDistribution(self, plot=True):
        """ Obtain and plot the equilibrium probabilities of each macrostate

        Parameters
        ----------
        plot : bool, optional, default=True
            Disable plotting of the probabilities by setting it to False

        Returns
        -------
        eq : ndarray
            An array of equilibrium probabilities of the macrostates

        Examples
        --------
        >>> model = Model(data)
        >>> model.markovModel(100, 5)
        >>> model.eqDistribution()
        """
        self._integrityCheck(postmsm=True)
        macroeq = np.ones(self.macronum) * -1
        for i in range(self.macronum):
            macroeq[i] = np.sum(self.msm.stationary_distribution[self.macro_ofmicro == i])

        if plot:
            from matplotlib import pylab as plt
            plt.ion()
            plt.figure()
            plt.bar(range(self.macronum), macroeq)
            plt.ylabel('Equilibrium probability')
            plt.xlabel('Macrostates')
            plt.xticks(np.arange(0.4, self.macronum+0.4, 1), range(self.macronum))
            plt.show()
        return macroeq

def plot_fullstack( binning = np.linspace(0,10,1), myquery='', plotvar = default_plot_variable, \
                    scalefactor = 1., user_ylim = None):

    fig = plt.figure(figsize=(10,6))
    plt.grid(True)
    lasthist = 0
    myhistos = gen_histos(binning=binning,myquery=myquery,plotvar=plotvar,scalefactor=scalefactor)
    for key, (hist, bins) in myhistos.iteritems():

      plt.bar(bins[:-1],hist,
              width=bins[1]-bins[0],
              color=colors[key],
              bottom = lasthist,
              edgecolor = 'k',
              label='%s: %d Events'%(labels[key],sum(hist)))
      lasthist += hist
     

    plt.title('CCSingleE Stacked Backgrounds',fontsize=25)
    plt.ylabel('Events',fontsize=20)
    if plotvar == '_e_nuReco' or plotvar == '_e_nuReco_better':
        xstring = 'Reconstructed Neutrino Energy [GeV]' 
    elif plotvar == '_e_CCQE':
        xstring = 'CCQE Energy [GeV]'
    else:
        xstring = plotvar
    plt.xlabel(xstring,fontsize=20)
    plt.legend()
    plt.xticks(list(plt.xticks()[0]) + [binning[0]])
    plt.xlim([binning[0],binning[-1]])

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 main(argv=None):
    if argv is None:
        argv = sys.argv
    if len(argv) != 4:
        print "Usage: " + argv[0] + " <report name> <report dir> <graph dir>"
        return 1

    report = argv[1]
    report_dir = argv[2]
    graph_dir = argv[3]

    file_list = glob(os.path.join(report_dir, 'part*'))
    #Copy the raw data file to the graph_dir
    raw_file = os.path.join(graph_dir, report + '.tsv')
    shutil.copyfile(file_list[0], raw_file)

    #Process the file into a graph, ideally I would combine the two into one but for now I'll stick with two
    data_file = csv.DictReader(open(raw_file, 'rb'), fieldnames = ['hour', 'Requests', 'Bytes'], delimiter="\t")
    
    #Make an empty set will all the hours in it os if an hour is not in the data it will be 0
    length = 24
    requests_dict = {}
    for num in range(length):
        requests_dict['%0*d' % (2, num)] = (0, 0)

    #add the values we have to the dictionaries
    for row in data_file:
        requests_dict[row['hour']] = (int(row['Requests']), int(row['Bytes']))

    #Now get the lists for graphing in the right order
    requests = []
    num_bytes = []
    requests_lists = requests_dict.items()
    requests_lists.sort(key=lambda req: req[0])
    for req in requests_lists:
        requests.append(req[1][0])
        num_bytes.append(req[1][1])

    fig = pylab.figure(1)
    pos = pylab.arange(length) + .5
    pylab.bar(pos, requests[:length], aa=True, ecolor='r')
    pylab.ylabel('Requests')
    pylab.xlabel('Hour')
    pylab.title('Request per hour')
    pylab.grid(True)

    #Save the figure
    pylab.savefig(os.path.join(graph_dir, report + '_requests.pdf'), bbox_inches='tight', pad_inches=1)

    #bytes listed 
    fig = pylab.figure(2)
    pos = pylab.arange(length) + .5
    pylab.bar(pos, num_bytes[:length], aa=True, ecolor='r')
    pylab.ylabel('Bytes')
    pylab.xlabel('Hour')
    pylab.title('Bytes per hour')
    pylab.grid(True)

    #Save the figure
    pylab.savefig(os.path.join(graph_dir, report + '_bytes.pdf'), bbox_inches='tight', pad_inches=1)

    def plot_histogram(self, main="", numrows=1, numcols=1, fignum=1):
        """Plot a histogram of choices and probability sums. Expects probabilities as (at least) a 2D array.
        """
        from matplotlib.pylab import bar, xticks, yticks, title, text, axis, figure, subplot

        probabilities = self.get_probabilities()
        if probabilities.ndim < 2:
            raise StandardError, "probabilities must have at least 2 dimensions."
        alts = probabilities.shape[1]
        width_par = (1 / alts + 1) / 2.0
        choice_counts = self.get_choice_histogram(0, alts)
        sum_probs = self.get_probabilities_sum()

        subplot(numrows, numcols, fignum)
        bar(arange(alts), choice_counts, width=width_par)
        bar(arange(alts) + width_par, sum_probs, width=width_par, color="g")
        xticks(arange(alts))
        title(main)
        Axis = axis()
        text(
            alts + 0.5,
            -0.1,
            "\nchoices histogram (blue),\nprobabilities sum (green)",
            horizontalalignment="right",
            verticalalignment="top",
        )

    def show(self,H):
        """
        Display the histogram of the data, together with the mixture model

        Parameters
        ----------
        H : ndarray
            The histogram of the data.
        """
        xmax = np.size(H)
        sH = np.sum(H)
        nH = H.astype(float)/sH

        L = np.zeros(xmax)
        ndraw = xmax-1
        for i in range(xmax):
            L0 = np.exp(self._bcoef(ndraw,i)+i*np.log(self.r0)+ \
                 (ndraw-i)*np.log(1-self.r0))
            L1 = np.exp(self._bcoef(ndraw,i)+i*np.log(self.r1)+ \
                 (ndraw-i)*np.log(1-self.r1))
            L[i] = self.Lambda*L0 + (1-self.Lambda)*L1

        L = L/L.sum()
        import matplotlib.pylab as mp
        mp.figure()
        mp.bar(np.arange(xmax),nH)
        mp.plot(np.arange(xmax)+0.5,L,'k',linewidth=2)

def Importance_Plot(data,label=None):
    '''
    :param data: DATAFRAME style
    :param label: y vector
    :param threshold: jude threshold
    :return: figure
    '''
    import numpy as np
    import matplotlib.pylab as plt
    from sklearn.ensemble import ExtraTreesClassifier
    import pandas as pd
    model=ExtraTreesClassifier()
    data1=np.array(data)
    model.fit(data1,label)
    importance=model.feature_importances_
    std = np.std([importance for tree in model.estimators_],axis=0)
    indices = np.argsort(importance)[::-1]
    namedata=data
    # Print the feature ranking
    print("Feature ranking:")
    importa=pd.DataFrame({'importance':list(importance[indices]),'Feature name':list(namedata.columns[indices])})
    print importa
    # Plot the feature importances of the forest
    plt.figure(figsize=(20, 8))
    plt.title("Feature importances")
    plt.bar(range(data1.shape[1]), importance[indices],
            color="g", yerr=std[indices], align="center")
    plt.xticks(range(data1.shape[1]), indices)
    plt.xlim([-1, data1.shape[1]])
    plt.grid(True)
    plt.show()

def plot_tuning_curves(direction_rates, title):
    """
    This function takes the x-values and the y-values  in units of spikes/s 
    (found in the two columns of direction_rates) and plots a histogram and 
    polar representation of the tuning curve. It adds the given title.
    """
        # yank columns and keep in correspondance 
    directions = direction_rates[0:][:,0]
    rates = direction_rates[0:][:,1]
    # histogram plot 
    plt.subplot(2, 2, 1)
    plt.title('Histogram ' + title)
    plt.axis([0, 360, 0, 70])
    plt.xlim(-22.5,337.5)
    plt.xlabel('Directions (in degrees)')
    plt.ylabel('Average Firing Rate (in spikes/s)')
    plt.bar(directions, rates, width=45, align='center')
    plt.xticks(directions)
    # polar plot 
    plt.subplot(2,2,2,polar=True)
    plt.title('Polar plot ' + title)
    #plt.legend('Diring Rate (spikes/s)')
    rates = np.append(rates, rates[0]) 
    theta = np.arange(0,361,45)*np.pi/180
    plt.polar(theta, rates)

    plt.show()
    # end plot_tuning_curves
    return 0

def plotHousing(impression):
    """
    生成房价随时间变化的图标
    """
    f = open("midWestHousingPrices.txt", 'r')
    #文件每一行是年季度价格
    labels, prices = [], []
    for line in f:
        year, quarter, price = line.split()
        label = year[2:4] + "\n Q" + quarter[1]
        labels.append(label)
        prices.append(float(price)/1000)
    #柱的X坐标
    quarters = np.arange(len(labels))
    #柱宽
    width = 0.5
    if impression == 'flat':
        plt.semilogy()
    plt.bar(quarters, prices, width, color='r')
    plt.xticks(quarters + width / 2.0, labels)
    plt.title("美国中西部各州房价")
    plt.xlabel("季度")
    plt.ylabel("平均价格($1000)")

    if impression == 'flat':
        plt.ylim(10, 10**3)
    elif impression == "volatile":
        plt.ylim(180, 220)
    elif impression == "fair":
        plt.ylim(150, 250)
    else:
        raise ValueError("Invalid input.")

def distributionValues(y,bins=None):
    if bins==None:
        bins=int(sqrt(len(y)))
    
    hist,bins=np.histogram(y,bins=bins)

    plt.bar(bins[1:],hist,width=bins[1]-bins[0])

    def plot_histogram_with_capacity(self, capacity, main=""):
        """Plot histogram of choices and capacities. The number of alternatives is determined
        from the second dimension of probabilities.
        """
        from matplotlib.pylab import bar, xticks, yticks, title, text, axis, figure, subplot

        probabilities = self.get_probabilities()
        if probabilities.ndim < 2:
            raise StandardError, "probabilities must have at least 2 dimensions."
        alts = self.probabilities.shape[1]
        width_par = (1 / alts + 1) / 2.0
        choice_counts = self.get_choice_histogram(0, alts)
        sum_probs = self.get_probabilities_sum()

        subplot(212)
        bar(arange(alts), choice_counts, width=width_par)
        bar(arange(alts) + width_par, capacity, width=width_par, color="r")
        xticks(arange(alts))
        title(main)
        Axis = axis()
        text(
            alts + 0.5,
            -0.1,
            "\nchoices histogram (blue),\ncapacities (red)",
            horizontalalignment="right",
            verticalalignment="top",
        )

def plotDist(subplot, X, Y, label):
    pylab.grid()
    pylab.subplot(subplot)
    pylab.bar(X, Y, 0.05)
    pylab.ylabel(label)
    pylab.xticks(arange(len(X)), X)
    pylab.yticks(arange(0,1,0.1))

def test_probabilities(exp, n=1000):
    d = {}
    for i in range(n):
        foo = rp.parsex(exp)
        # foo = len(foo.replace(' ',''))
        if foo in d.keys():
            d[foo] += 1
        else:
            d[foo] = 1
    # lists = sorted(d.items())
    # x, y = zip(*lists)  # unpack a list of pairs into two tuples
    # print x
    # print y
    # for a, b in zip(x, y):
    #     plt.text(a,b, str("%s\n%s" % (a, b)))
    plt.xlabel("String length")
    plt.ylabel("Occurence")
    plt.title("Union: %s P=0.3 N=1000" % exp)
    # plt.plot(x, y)

    # For bar chart (use on Union)
    # See for labeling: https://stackoverflow.com/a/30229062
    l = sorted(d.items())
    x, y = zip(*l)
    plt.bar(range(len(y)), y, align="center")
    plt.xticks(range(len(x)), x)

    plt.show()

 def eqDistribution(self, plot=True):
     if plot:
         from matplotlib import pylab as plt
         plt.bar(range(self.macronum), self.hmm.stationary_distribution)
         #ax = plt.gca()
         #ax.set_xticklabels([str(a) for a in self.hmm.active_set])
         plt.xticks(np.arange(self.macronum)+0.4, [str(a) for a in self.hmm.active_set])
     return self.hmm.stationary_distribution

 def plot_height(self):
     """Plot the height of the non-leaves nodes
     """
     import matplotlib.pylab as mp
     mp.figure()
     sh = np.sort(self.height[self.isleaf() == False])
     n = np.sum(self.isleaf() == False)
     mp.bar(np.arange(n), sh)

def LE(G, dim, verbose=0, maxiter=1000):
    """
    Laplacian Embedding of the data
    returns the dim-dimensional LE of the graph G
    
    Parameters
    ----------
    G, nipy.neurospin.graph.WeightedGraph instance that represents the data
    dim=1, int, number of dimensions
    verbose=0, verbosity level
    maxiter=1000, int, maximum number of iterations of the algorithm 
    
    Returns
    -------
    chart, array of shape(G.V,dim)
    
    Note
    ----
    In fact the current implementation retruns
    what is now referred to a diffusion map at time t=1
    """
    n = G.V
    dim = np.minimum(dim,n)
    chart = nr.randn(G.V,dim+2)
    f = ff.Field(G.V,G.edges,G.weights,chart)
    LNorm,RNorm = local_sym_normalize(G)
    # note : normally Rnorm = Lnorm
    if verbose:
        print np.sqrt(np.sum((LNorm-RNorm)**2))/np.sum(LNorm)
    eps = 1.e-7

    f1 = np.zeros((G.V,dim+2))
    for i in range(maxiter):
        f.diffusion(10)
        f0 = Orthonormalize(f.field)
        f.field = f0
        if nl.norm(f1-f0)<eps:
            break
        else:
            f1 = f0

    if verbose:
        print i,nl.norm(f1-f0)
    f.diffusion()
    f0 = f.field
    
    U,S,V = nl.svd(f0,0)
    RNorm = np.reshape(np.sqrt(RNorm),(n,1))
    chart = S[:dim]*np.repeat(1./RNorm,dim,1)*U[:,1:dim+1]
    chart = chart/np.sqrt(np.sum(chart**2,0))
            
    if verbose:
        import matplotlib.pylab as mp
        mp.figure()
        mp.bar(np.arange(np.size(S)-1),np.sqrt(1-S[1:]))
        print "laplacian eigenvalues: ",1-S
        
    return chart

def plot_results(results, title, xlabels, ylabel="Success Rate"):
    '''Plot a bar graph of results'''
    ind = np.arange(len(results))
    width = 0.4
    plt.bar(ind, results, width, color="#1AADA4")
    plt.ylabel(ylabel)
    plt.ylim(ymax=100)
    plt.xticks(ind+width/2.0, xlabels)
    plt.title(title)