def check_models(self):
        plt.figure('Bandgap narrowing')
        Na = np.logspace(12, 20)
        Nd = 0.
        dn = 1e14
        temp = 300.

        for author in self.available_models():
            BGN = self.update(Na=Na, Nd=Nd, nxc=dn,
                              author=author,
                              temp=temp)

            if not np.all(BGN == 0):
                plt.plot(Na, BGN, label=author)

        test_file = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            'Si', 'check data', 'Bgn.csv')

        data = np.genfromtxt(test_file, delimiter=',', names=True)

        for name in data.dtype.names[1:]:
            plt.plot(
                data['N'], data[name], 'r--',
                label='PV-lighthouse\'s: ' + name)

        plt.semilogx()
        plt.xlabel('Doping (cm$^{-3}$)')
        plt.ylabel('Bandgap narrowing (K)')

        plt.legend(loc=0)

def make_corr1d_fig(dosave=False):
    corr = make_corr_both_hemi()
    lw=2; fs=16
    pl.figure(1)#, figsize=(8, 7))
    pl.clf()
    pl.xlim(4,300)
    pl.ylim(-400,+500)    
    lambda_titles = [r'$20 < \lambda < 30$',
                     r'$30 < \lambda < 40$',
                     r'$\lambda > 40$']
    colors = ['blue','green','red']
    for i in range(3):
        corr1d, rcen = corr_1d_from_2d(corr[i])
        ipdb.set_trace()
        pl.semilogx(rcen, corr1d*rcen**2, lw=lw, color=colors[i])
        #pl.semilogx(rcen, corr1d*rcen**2, 'o', lw=lw, color=colors[i])
    pl.xlabel(r'$s (Mpc)$',fontsize=fs)
    pl.ylabel(r'$s^2 \xi_0(s)$', fontsize=fs)    
    pl.legend(lambda_titles, 'lower left', fontsize=fs+3)
    pl.plot([.1,10000],[0,0],'k--')
    s_bao = 149.28
    pl.plot([s_bao, s_bao],[-9e9,+9e9],'k--')
    pl.text(s_bao*1.03, 420, 'BAO scale')
    pl.text(s_bao*1.03, 370, '%0.1f Mpc'%s_bao)
    if dosave: pl.savefig('xi1d_3bin.pdf')

 def plot_bernoulli_matrix(self, show_npfs=False):
   """
   Plot the heatmap of the Bernoulli matrix 
   @self
   @show_npfs - Highlight NPFS detections [Boolean] 
   """
   matrix = self.Bernoulli_matrix
   if show_npfs == False:
     plot = plt.imshow(matrix)
     plot.set_cmap('hot')
     plt.colorbar()
     plt.xlabel("Bootstraps")
     plt.ylabel("Feature")
     plt.show()
   else:
     for i in self.selected_features:
       for k in range(len(matrix[i])):
         matrix[i,k] = .5
     plot = plt.imshow(matrix)
     plot.set_cmap('hot')
     plt.xlabel("Bootstraps")
     plt.ylabel("Feature")
     plt.colorbar()
     plt.show()
   return None

    def bar(self, key_word_sep = " ", title=None, **kwargs):
        """Generates a pylab bar plot from the result set.

        ``matplotlib`` must be installed, and in an
        IPython Notebook, inlining must be on::

            %%matplotlib inline

        The last quantitative column is taken as the Y values;
        all other columns are combined to label the X axis.

        Parameters
        ----------
        title: Plot title, defaults to names of Y value columns
        key_word_sep: string used to separate column values
                      from each other in labels

        Any additional keyword arguments will be passsed
        through to ``matplotlib.pylab.bar``.
        """
        import matplotlib.pylab as plt
        self.guess_pie_columns(xlabel_sep=key_word_sep)
        plot = plt.bar(range(len(self.ys[0])), self.ys[0], **kwargs)
        if self.xlabels:
            plt.xticks(range(len(self.xlabels)), self.xlabels,
                       rotation=45)
        plt.xlabel(self.xlabel)
        plt.ylabel(self.ys[0].name)
        return plot

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_values(self, TITLE, SAVE):
     plot(self.list_of_densities, self.list_of_pressures)
     title(TITLE)
     xlabel("Densities")
     ylabel("Pressure")
     savefig(SAVE)
     show()

    def check_models(self):
        '''
        Displays a plot of the models against that taken from a
        respected website (https://www.pvlighthouse.com.au/)
        '''
        plt.figure('Intrinsic bandgap')
        t = np.linspace(1, 500)

        for author in self.available_models():

            Eg = self.update(temp=t, author=author, multiplier=1.0)
            plt.plot(t, Eg, label=author)

        test_file = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            'Si', 'check data', 'iBg.csv')

        data = np.genfromtxt(test_file, delimiter=',', names=True)

        for temp, name in zip(data.dtype.names[0::2], data.dtype.names[1::2]):
            plt.plot(
                data[temp], data[name], '--', label=name)

        plt.xlabel('Temperature (K)')
        plt.ylabel('Intrinsic Bandgap (eV)')

        plt.legend(loc=0)
        self.update(temp=0, author=author, multiplier=1.01)

def fdr(p_values=None, verbose=0):
    """Returns the FDR associated with each p value

    Parameters
    -----------
    p_values : ndarray of shape (n)
        The samples p-value

    Returns
    -------
    q : array of shape(n)
        The corresponding fdr values
    """
    p_values = check_p_values(p_values)
    n_samples = p_values.size
    order = p_values.argsort()
    sp_values = p_values[order]

    # compute q while in ascending order
    q = np.minimum(1, n_samples * sp_values / np.arange(1, n_samples + 1))
    for i in range(n_samples - 1, 0, - 1):
        q[i - 1] = min(q[i], q[i - 1])

    # reorder the results
    inverse_order = np.arange(n_samples)
    inverse_order[order] = np.arange(n_samples)
    q = q[inverse_order]

    if verbose:
        import matplotlib.pylab as mp
        mp.figure()
        mp.xlabel('Input p-value')
        mp.plot(p_values, q, '.')
        mp.ylabel('Associated fdr')
    return q

    def plot_corner_posteriors(self, savefile=None, labels=["T1", "R1", "Av", "T2", "R2"]):
        '''
        Plots the corner plot of the MCMC results.
        '''
        ndim = len(self.sampler.flatchain[0,:])
        chain = self.sampler
        samples = chain.flatchain
        
        samples = samples[:,0:ndim]  
        plt.figure(figsize=(8,8))
        fig = corner.corner(samples, labels=labels[0:ndim])
        plt.title("MJD: %.2f"%self.mjd)
        name = self._get_save_path(savefile, "mcmc_posteriors")
        plt.savefig(name)
        plt.close("all")
        

        plt.figure(figsize=(8,ndim*3))
        for n in range(ndim):
            plt.subplot(ndim,1,n+1)
            chain = self.sampler.chain[:,:,n]
            nwalk, nit = chain.shape
            
            for i in np.arange(nwalk):
                plt.plot(chain[i], lw=0.1)
                plt.ylabel(labels[n])
                plt.xlabel("Iteration")
        name_walkers = self._get_save_path(savefile, "mcmc_walkers")
        plt.tight_layout()
        plt.savefig(name_walkers)
        plt.close("all")  

def plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):
    """
    Plot a radiallysymmetric Q model.

    plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):

    r_min=minimum radius [km], r_max=maximum radius [km], dr=radius
    increment [km]

    Currently available models (model): cem, prem, ql6
    """
    import matplotlib.pylab as plt

    r = np.arange(r_min, r_max + dr, dr)
    q = np.zeros(len(r))

    for k in range(len(r)):

        if model == 'cem':
            q[k] = q_cem(r[k])
        elif model == 'ql6':
            q[k] = q_ql6(r[k])
        elif model == 'prem':
            q[k] = q_prem(r[k])

    plt.plot(r, q, 'k')
    plt.xlim((0.0, r_max))
    plt.xlabel('radius [km]')
    plt.ylabel('Q')
    plt.show()

    def handle(self, *args, **options):
        try:
            from matplotlib import pylab as pl
            import numpy as np
        except ImportError:
            raise Exception('Be sure to install requirements_scipy.txt before running this.')

        all_names_and_counts = RawCommitteeTransactions.objects.all().values('attest_by_name').annotate(total=Count('attest_by_name')).order_by('-total')
        all_names_and_counts_as_tuple_and_sorted = sorted([(row['attest_by_name'], row['total']) for row in all_names_and_counts], key=lambda row: row[1])
        print "top ten attestors:  (name, number of transactions they attest for)"
        for row in all_names_and_counts_as_tuple_and_sorted[-10:]:
            print row

        n_bins = 100
        filename = 'attestor_participation_distribution.png'

        x_max = all_names_and_counts_as_tuple_and_sorted[-31][1]  # eliminate top outliers from hist
        x_min = all_names_and_counts_as_tuple_and_sorted[0][1]

        counts = [row['total'] for row in all_names_and_counts]
        pl.figure(1, figsize=(18, 6))
        pl.hist(counts, bins=np.arange(x_min, x_max, (float(x_max)-x_min)/100) )
        pl.title('Histogram of Attestor Participation in RawCommitteeTransactions')
        pl.xlabel('Number of transactions a person attested for')
        pl.ylabel('Number of people')
        pl.savefig(filename)

def flipPlot(minExp, maxExp):
    """假定minEXPy和maxExp是正整数且minExp<maxExp
    绘制出2**minExp到2**maxExp次抛硬币的结果
    """
    ratios = []
    diffs = []
    aAxis = []
    for i in range(minExp, maxExp+1):
        aAxis.append(2**i)
    for numFlips in aAxis:
        numHeads = 0
        for n in range(numFlips):
            if random.random() < 0.5:
                numHeads += 1
        numTails = numFlips - numHeads
        ratios.append(numHeads/numFlips)
        diffs.append(abs(numHeads-numTails))
    plt.figure()
    ax1 = plt.subplot(121)
    plt.title("Difference Between Heads and Tails")
    plt.xlabel('Number of Flips')
    plt.ylabel('Abs(#Heads - #Tails)')
    ax1.semilogx(aAxis, diffs, 'bo')
    ax2 = plt.subplot(122)
    plt.title("Heads/Tails Ratios")
    plt.xlabel('Number of Flips')
    plt.ylabel("#Heads/#Tails")
    ax2.semilogx(aAxis, ratios, 'bo')
    plt.show()

def plot_runtime_results(results):
    plt.rcParams["figure.figsize"] = 7,7
    plt.rcParams["font.size"] = 22
    matplotlib.rc("xtick", labelsize=24)
    matplotlib.rc("ytick", labelsize=24)

    params = {"text.fontsize" : 32,
              "font.size" : 32,
              "legend.fontsize" : 30,
              "axes.labelsize" : 32,
              "text.usetex" : False
              }
    plt.rcParams.update(params)
    
    #plt.semilogx(results[:,0], results[:,3], 'r-x', lw=3)
    #plt.semilogx(results[:,0], results[:,1], 'g-D', lw=3)
    #plt.semilogx(results[:,0], results[:,2], 'b-s', lw=3)

    plt.plot(results[:,0], results[:,3], 'r-x', lw=3, ms=10)
    plt.plot(results[:,0], results[:,1], 'g-D', lw=3, ms=10)
    plt.plot(results[:,0], results[:,2], 'b-s', lw=3, ms=10)

    plt.legend(["Chain", "Tree", "FFT Tree"], loc="upper left")
    plt.xticks([1e5, 2e5, 3e5])
    plt.yticks([0, 60, 120, 180])

    plt.xlabel("Problem Size")
    plt.ylabel("Runtime (sec)")
    return results

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 inter_show(start, lc, eta, vol_ins, props, lbl_outs, grdts, pars):
    '''
    Plots a display of training information to the screen
    '''
    import matplotlib.pylab as plt
    name_in, vol  = vol_ins.popitem()
    name_p,  prop = props.popitem()
    name_l,  lbl  = lbl_outs.popitem()
    name_g,  grdt = grdts.popitem()

    m_input = volume_util.crop(vol[0,:,:,:], prop.shape[-3:]) #good enough for now

    # real time visualization
    plt.subplot(251),   plt.imshow(vol[0,0,:,:],    interpolation='nearest', cmap='gray')
    plt.xlabel('input')
    plt.subplot(252),   plt.imshow(m_input[0,:,:],    interpolation='nearest', cmap='gray')
    plt.xlabel('matched input')
    plt.subplot(253),   plt.imshow(prop[0,0,:,:],   interpolation='nearest', cmap='gray')
    plt.xlabel('output')
    plt.subplot(254),   plt.imshow(lbl[0,0,:,:],    interpolation='nearest', cmap='gray')
    plt.xlabel('label')
    plt.subplot(255),   plt.imshow(grdt[0,0,:,:],   interpolation='nearest', cmap='gray')
    plt.xlabel('gradient')

    plt.subplot(256)
    plt.plot(lc.tn_it, lc.tn_err, 'b', label='train')
    plt.plot(lc.tt_it, lc.tt_err, 'r', label='test')
    plt.xlabel('iteration'), plt.ylabel('cost energy')
    plt.subplot(257)
    plt.plot( lc.tn_it, lc.tn_cls, 'b', lc.tt_it, lc.tt_cls, 'r')
    plt.xlabel('iteration'), plt.ylabel( 'classification error' )
    return

def modelfit(alg, train, target, test, useTrainCV=True, cv_folds=5, early_stopping_rounds=50):
    
    if useTrainCV:
        xgboost_params = alg.get_xgb_params()
        xgtrain = xgb.DMatrix(train.values, label=target.values)
        xgtest = xgb.DMatrix(test.values)
        watchlist = [(xgtrain, 'train')] # Specify validations set to watch performance
        cvresult = xgb.cv(xgboost_params, xgtrain, num_boost_round=alg.get_params()['n_estimators'], nfold=cv_folds, early_stopping_rounds=early_stopping_rounds) #metrics='auc',show_progress=False
        alg.set_params(n_estimators=cvresult.shape[0])
    
    # Fit the algorithm on the data
    alg.fit(train, target, eval_metric='auc')

    # Predict training set:
    train_preds = alg.predict(train)
    train_predprob = alg.predict_proba(train)[:,1]
    
    # Print model report:
    print "\nModel Report"
    print "Accuracy : %.4g" % metrics.accuracy_score(target.values, train_preds)
    print "AUC Score (Train): %f" % metrics.roc_auc_score(target, train_predprob)

    # Make a prediction:
    print('Predicting......')
    test_predprob = alg.predict_proba(test)[:,1]

    feat_imp = pd.Series(alg.booster().get_fscore()).sort_values(ascending=False)
    feat_imp.plot(kind='bar', title='Feature Importances')
    plt.ylabel('Feature Importance Score')
    #plt.show()
    return test_predprob

def plot_confusion_matrix(cm, title='', cmap=plt.cm.Blues):
    #print cm
    #display vehicle, idle, walking accuracy respectively
    #display overall accuracy
    print type(cm)
   # plt.figure(index
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    #plt.figure("")
    plt.title("Confusion Matrix")
    plt.colorbar()
    tick_marks = [0,1,2]
    target_name = ["driving","idling","walking"]


    plt.xticks(tick_marks,target_name,rotation=45)

    plt.yticks(tick_marks,target_name,rotation=45)
    print len(cm[0])

    for i in range(0,3):
        for j in range(0,3):
         plt.text(i,j,str(cm[i,j]))
    plt.tight_layout()
    plt.ylabel("Actual Value")
    plt.xlabel("Predicted Outcome")

	def plot_cell(self,cell_number=0,label='insert_label'):

		current_cell = self.cell_list[cell_number]
		temp = current_cell.temp
		cd_signal = current_cell.cd_signal
		cd_calc = current_cell.cd_calc()
		
		ax = pylab.gca()

		pylab.plot(temp,cd_signal,'o',color='black')
                pylab.plot(temp,cd_calc,color='black')
		pylab.xlabel(r'Temperature ($^{\circ}$C)')
		pylab.ylabel('mdeg')
		pylab.ylim([-25,-4])
		dH = numpy.round(current_cell.dH, decimals=1)
		Tm = numpy.round(current_cell.Tm-273.15, decimals=1)
		nf = current_cell.nf
		nu = current_cell.nu
		textstr_dH = '${\Delta}H_{m}$ = %.1f kcal/mol' %dH
		textstr_Tm ='$T_{m}$ = %.1f $^{\circ}$C' %Tm
		textstr_nf ='$N_{folded}$ = %d' %nf
		textstr_nu ='$N_{unfolded}$ = %d'%nu
		ax.text(8,-6,textstr_dH, fontsize=16,ha='left',va='top')
		ax.text(8,-7.5,textstr_Tm, fontsize=16,ha='left',va='top')
		ax.text(8,-9,textstr_nf, fontsize=16,ha='left',va='top')
		ax.text(8,-10.5,textstr_nu, fontsize=16,ha='left',va='top')
		pylab.title(label)		
		pylab.show()

		return

 def test_simple_gen(self):
     self_con = .8
     other_con = 0.05
     g = self.gen.gen_stoch_blockmodel(min_degree=1, blocks=5, self_con=self_con, other_con=other_con,
                                       powerlaw_exp=2.1, degree_seq='powerlaw', num_nodes=1000, num_links=3000)
     deg_hist = vertex_hist(g, 'total')
     res = fit_powerlaw.Fit(g.degree_property_map('total').a, discrete=True)
     print 'powerlaw alpha:', res.power_law.alpha
     print 'powerlaw xmin:', res.power_law.xmin
     if len(deg_hist[0]) != len(deg_hist[1]):
         deg_hist[1] = deg_hist[1][:len(deg_hist[0])]
     print 'plot degree dist'
     plt.plot(deg_hist[1], deg_hist[0])
     plt.xscale('log')
     plt.xlabel('degree')
     plt.ylabel('#nodes')
     plt.yscale('log')
     plt.savefig('deg_dist_test.png')
     plt.close('all')
     print 'plot graph'
     pos = sfdp_layout(g, groups=g.vp['com'], mu=3)
     graph_draw(g, pos=pos, output='graph.png', output_size=(800, 800),
                vertex_size=prop_to_size(g.degree_property_map('total'), mi=2, ma=30), vertex_color=[0., 0., 0., 1.],
                vertex_fill_color=g.vp['com'],
                bg_color=[1., 1., 1., 1.])
     plt.close('all')
     print 'init:', self_con / (self_con + other_con), other_con / (self_con + other_con)
     print 'real:', gt_tools.get_graph_com_connectivity(g, 'com')

def plot_peaks(spectra, ref, zl, pl, filename=''):
  plt.figure();
  plt.title("Debug: Peak fitting for '%s'" % filename);
  plt.xlabel("y-position [px]");
  plt.ylabel("Intensity");

  Nspectra, Npx = spectra.shape;

  for s in xrange(Nspectra):
    scale = 1./spectra.max();
    offset= -s*0.1;

    # plot data
    plt.plot(spectra[s]*scale + offset,'k',linewidth=2);

    # plot first peak
    p,A,w = zl[s];
    x = np.arange(-2*w, 2*w) + p;
    plt.plot(x,gauss(x,*zl[s])*scale + offset,'r');
  
    # plot second peak
    if ref is not None:
      p,A   = pl[s];
      x = np.arange(len(ref)) - len(ref)/2 + p;
      plt.plot(x,ref/ref.max()*A*scale + offset,'g');