def _create_histogram(M_c, data, columns, mc_col_indices, filename):
    dir = S.path.web_resources_data_dir
    full_filename = os.path.join(dir, filename)
    num_rows = data.shape[0]
    num_cols = data.shape[1]

    p.figure()
    # col_i goes from 0 to number of predicted columns
    # mc_col_idx is the original column's index in M_c
    for col_i in range(num_cols):
        mc_col_idx = mc_col_indices[col_i]
        data_i = data[:, col_i]
        ax = p.subplot(1, num_cols, col_i, title=columns[col_i])
        if M_c['column_metadata'][mc_col_idx]['modeltype'] == 'normal_inverse_gamma':
            p.hist(data_i, orientation='horizontal')
        else:
            str_data = [du.convert_code_to_value(M_c, mc_col_idx, code) for code in data_i]
            unique_labels = list(set(str_data))
            np_str_data = np.array(str_data)
            counts = []
            for label in unique_labels:
                counts.append(sum(np_str_data == label))
            num_vals = len(M_c['column_metadata'][mc_col_idx]['code_to_value'])
            rects = p.barh(range(num_vals), counts)
            heights = np.array([rect.get_height() for rect in rects])
            ax.set_yticks(np.arange(num_vals) + heights/2)
            ax.set_yticklabels(unique_labels)

    p.tight_layout()
    p.savefig(full_filename)

    def plot_number_alteration_by_tissue(self, fontsize=10, width=0.9):
        """Plot number of alterations

        .. plot::
            :width: 100%
            :include-source:

            from gdsctools import *
            data = gdsctools_data("test_omnibem_genomic_alterations.csv.gz")
            bem = OmniBEMBuilder(data)
            bem.filter_by_gene_list(gdsctools_data("test_omnibem_genes.txt"))
            bem.plot_number_alteration_by_tissue()

        """
        count = self.unified.groupby(['TISSUE_TYPE'])['GENE'].count()
        try:
            count.sort_values(inplace=True, ascending=False)
        except:
            count.sort(inplace=True, ascending=False)
        count.plot(kind="bar", width=width)
        pylab.grid()
        pylab.xlabel("Tissue Type", fontsize=fontsize)
        pylab.ylabel("Total number of alterations in cell lines",
                     fontsize=fontsize)
        try:pylab.tight_layout()
        except:pass

def plot_Barycenter(dataset_name, feat, unfeat, repo):

    if dataset_name==MNIST:
        _, _, test=get_data(dataset_name, repo, labels=True)
        xtest1,_,_, labels,_=test
    else:
        _, _, test=get_data(dataset_name, repo, labels=False)
        xtest1,_,_ =test
        labels=np.zeros((len(xtest1),))
    # get labels
    def bary_wdl2(index): return _bary_wdl2(index, xtest1, feat, unfeat)
    
    n=xtest1.shape[-1]
    
    num_class = (int)(max(labels)+1)
    barys=[bary_wdl2(np.where(labels==i)) for i in range(num_class)]
    pl.figure(1, (num_class, 1))
    for i in range(num_class):
        pl.subplot(1,10,1+i)
        pl.imshow(barys[i][0,0,:,:],cmap='Blues',interpolation='nearest')
        pl.xticks(())
        pl.yticks(())
        if i==0:
            pl.ylabel('DWE Bary.')
        if num_class >1:
            pl.title('{}'.format(i))
    pl.tight_layout(pad=0,h_pad=-2,w_pad=-2) 
    pl.savefig("imgs/{}_dwe_bary.pdf".format(dataset_name))

def param_sweeping(clf, obj, X, y, param_dist, metric, param, clfName):
	'''Plot a parameter sweeping (ps) curve with the param_dist as a axis, and the scoring based on metric as y axis.
	Keyword arguments:
	clf - - classifier
	X - - feature matrix
	y - - target array
	param - - a parameter of the classifier
	param_dist - - the parameter distribution of param
	clfName - - the name of the classifier
	metric - - the metric we use to evaluate the performance of the classifiers
	obj - - the name of the dataset we are using'''
	scores = []
	for i in param_dist:
		y_true = []
		y_pred = []
		# new classifer each iteration
		newclf = eval("clf.set_params("+ param + "= i)")
		y_pred, y_true, gs_score_list, amp = testAlgo(newclf, X, y, clfName)
		mean_fpr, mean_tpr, mean_auc = plot_unit_prep(y_pred, y_true, metric)
		scores.append(mean_auc)
		print("Area under the ROC curve : %f" % mean_auc)
	fig = pl.figure(figsize=(8,6),dpi=150)
	paramdist_len = len(param_dist)
	pl.plot(range(paramdist_len), scores, label = 'Parameter Sweeping Curve')
	pl.xticks(range(paramdist_len), param_dist, size = 15, rotation = 45)
	pl.xlabel(param.upper(),fontsize=30)
	pl.ylabel(metric.upper(),fontsize=30)
	pl.title('Parameter Sweeping Curve',fontsize=25)
	pl.legend(loc='lower right')
	pl.tight_layout()
	fig.savefig('plots/'+obj+'/'+ clfName +'_' + param +'_'+'ps.pdf')
	pl.show()

def plotMixNB(X, y, obj, featureNames, whichMix):
	"""Plot MixNB fit on top of X.
	"""
	save_path = '../MSPrediction-Python/plots/'+obj+'/'+whichMix
	clf = classifiers[whichMix]
	clf,_,_ = fitAlgo(clf, X,y, opt= True, param_dict = param_dist_dict[whichMix])
	unique_y = np.unique(y)
	# norm_func = lambda x, sigma, theta: 1 if np.isnan(x) else -0.5 * np.log(2 * np.pi*sigma) - 0.5 * ((x - theta)**2/sigma) 
	# norm_func = np.vectorize(norm_func)
	n_samples = X.shape[0]
	for j in range(X.shape[1]):
		fcol = X[:,j]
		optmodel = clf.optmodels[:,j]
		distname = clf.distnames[j]
		jfeature = featureNames[j]
		jpath = save_path +'_'+jfeature+'.pdf'
		fig = pl.figure(figsize=(8,6),dpi=150)
		for i, y_i in enumerate(unique_y):
			fcoli = fcol[y == y_i]
			itfreq = itemfreq(fcoli)
			uniqueVars = itfreq[:,0]
			freq = itfreq[:,1]
			freq = freq/sum(freq)
			pred = np.exp(optmodel[i](uniqueVars))
			# print pred
			# print pred
			pl.plot(uniqueVars, pred, label= str(y_i)+'_model')
			pl.plot(uniqueVars, freq, label= str(y_i) +'_true')
		pl.xlabel(jfeature)
		pl.ylabel("density")
		pl.title(distname)
		pl.legend(loc='best')
		pl.tight_layout()
		# pl.show()
		fig.savefig(jpath)

def pca(inF,MIN):
    df = pd.read_table(inF, header=0)
    dc = list(df.columns)
    dc[0]='GeneID'
    df.columns = dc
    
    print(df.shape)
    sel = True 
    for i in range(4, df.shape[1]-1):
        sel = (df.ix[:,i] < MIN) & (df.ix[:,i+1]< MIN)
    df = df.ix[~sel,:]
    print(df.shape)
    
    X = df.ix[:,4:df.shape[1]].values.T
    y = df.columns[4:df.shape[1]].values
    X_std = StandardScaler().fit_transform(X)
    
    #pca = PCA(n_components=2)
    pca = PCA()
    Y_sklearn = pca.fit_transform(X_std)
    
    
    fig = plt.figure()
    plt.style.use('ggplot')
    #plt.style.use('seaborn-whitegrid')
    ax = fig.add_subplot(111)
    for lab, col in zip(y,['r','g','b','c'] + sns.color_palette("cubehelix", df.shape[1]-4-4)):
        ax.scatter(Y_sklearn[y==lab, 0],Y_sklearn[y==lab, 1],label=lab,c=col, s=80)
    
    
    ax.set_xlabel('Principal Component 1 : %.2f'%(pca.explained_variance_ratio_[0]*100) + '%')
    ax.set_ylabel('Principal Component 2 : %.2f'%(pca.explained_variance_ratio_[1]*100) + '%')
    ax.legend(loc='upper right', prop={'size':8})
    plt.tight_layout()
    plt.savefig(inF + '-MIN' + str(MIN) + '.pdf')

def tracks_movie(base, skip=1, frames=500, size=10):
    """
    A movie of each particle as a point
    """
    conf, track, pegs = load(base)

    fig = pl.figure(figsize=(size,size*conf['top']/conf['wall']))
    plot = None

    for t in xrange(1,max(frames, track.shape[1]/skip)):
        tmp = track[:,t*skip,:]
        if not ((tmp[:,0] > 0) & (tmp[:,1] > 0) & (tmp[:,0] < conf['wall']) & (tmp[:,1] < conf['top'])).any():
            continue

        if plot is None:
            plot = pl.plot(tmp[:,0], tmp[:,1], 'k,', alpha=1.0, ms=0.1)[0]
            pl.xticks([])
            pl.yticks([])
            pl.xlim(0,conf['wall'])
            pl.ylim(0,conf['top'])
            pl.tight_layout()
        else:
            plot.set_xdata(tmp[:,0])
            plot.set_ydata(tmp[:,1])
        pl.draw()
        pl.savefig(base+'-movie-%05d.png' % (t-1))

  def make_plotII(self):
    # retrieve data
    D=self.D

    kmap={}
    kmap['Q2 = 2']   = {'c':'r','ls':'-'}
    kmap['Q2 = 5']   = {'c':'g','ls':'--'}
    kmap['Q2 = 10']  = {'c':'b','ls':'-.'}
    kmap['Q2 = 100'] = {'c':'k','ls':':'}


    ax=py.subplot(111)
    DF=D['AV18']
    for Q2 in [2,5,10,100]:
      k='Q2 = %d'%Q2
      Q2=float(k.split('=')[1])
      DF=D['AV18'][D['AV18'].Q2==Q2]
      cls=kmap[k]['c']+kmap[k]['ls']
      ax.plot(DF.X,DF.THEORY,cls,lw=2.0,label=r'$Q^2=%0.0f~{\rm GeV}^2$'%Q2)

    ax.set_xlabel('$x$',size=25)
    ax.set_ylabel(r'$F_2^d\, /\, F_2^N$',size=25)
    ax.set_ylim(0.97,1.08)
    ax.axhline(1,color='k',ls='-',alpha=0.2)

    ax.legend(frameon=0,loc=2,fontsize=22)
    py.tick_params(axis='both',labelsize=22)
    py.tight_layout()
    py.savefig('gallery/F2d_F2_II.pdf')

     def run(self):
         plts = {}
         graphs = {}
         pos  = 0

         plt.ion()
         plt.style.use('ggplot')

         for name in sorted(self.names.values()):
             p = plt.subplot(math.ceil(len(self.names) / 2), 2, pos+1)
             p.set_ylim([0, 100])
             p.set_title(self.machine_classes[name] + " " + name)

             p.get_xaxis().set_visible(False)

             X = range(0, NUM_ENTRIES, 1)
             Y = NUM_ENTRIES * [0]
             graphs[name] = p.plot(X, Y)[0]

             plts[name] = p
             pos += 1 

         plt.tight_layout()
         
         while True:
             for name, p in plts.items():
                 graphs[name].set_ydata(self.loads[name])

             plt.draw()
             plt.pause(0.05)

def runSimulation1(numSteps):
    """
    Runs the simulation for `numSteps` time steps.

    Returns a tuple of two lists: (rabbit_populations, fox_populations)
      where rabbit_populations is a record of the rabbit population at the
      END of each time step, and fox_populations is a record of the fox population
      at the END of each time step.

    Both lists should be `numSteps` items long.
    """

    rabbit_sim = []
    fox_sim = []
    for step in range(numSteps):
        rabbitGrowth()
        rabbit_sim.append(CURRENTRABBITPOP)
        foxGrowth()
        fox_sim.append(CURRENTFOXPOP)
    #print "CURRENTRABBITPOP", CURRENTRABBITPOP
    #print "CURRENTFOXPOP", CURRENTFOXPOP
    pylab.plot(range(numSteps), rabbit_sim, '-g', label='Rabbit population')
    pylab.plot(range(numSteps), fox_sim, '-o', label='Fox population')
    pylab.title('Fox and rabbit population in the wood')
    xlabel = "Plot for simulation of {} steps".format(numSteps)
    pylab.xlabel(xlabel)
    pylab.ylabel('Current fox and rabbit population')
    pylab.legend(loc='upper right')
    pylab.tight_layout()
    pylab.show()
    pylab.clf()

  def make_plotI(self):
    # retrieve data
    D=self.D

    kmap={}
    kmap['AV18']   = {'c':'r','ls':'-'}
    kmap['CDBONN'] = {'c':'g','ls':'--'}
    kmap['WJC1']   = {'c':'k','ls':'-.'}
    kmap['WJC2']   = {'c':'b','ls':':'}

    ax=py.subplot(111)
    for k in ['AV18','CDBONN','WJC1','WJC2']:
      DF=D[k]
      DF=DF[DF.Q2==10]
      if k=='CDBONN':
        label='CDBonn'
      else:
        label=k
      cls=kmap[k]['c']+kmap[k]['ls']
      ax.plot(DF.X,DF.THEORY,cls,lw=2.0,label=tex(label))

    ax.set_xlabel('$x$',size=25)
    ax.set_ylabel(r'$F_2^d\, /\, F_2^N$',size=25)
    ax.set_ylim(0.97,1.08)
    ax.axhline(1,color='k',ls='-',alpha=0.2)

    ax.legend(frameon=0,loc=2,fontsize=22)
    py.tick_params(axis='both',labelsize=22)
    py.tight_layout()
    py.savefig('gallery/F2d_F2_I.pdf')
    py.close()

def createMP4(text, filename, thumbnail, faimsSync):
	fig = pylab.plt.figure()
	fig.suptitle(text)
	ax = fig.add_subplot(111)
	ax.set_aspect('equal')
	ax.get_xaxis().set_visible(False)
	ax.get_yaxis().set_visible(False)

	im = ax.imshow(numpy.random.rand(300,300),cmap='gray',interpolation='nearest')
	im.set_clim([0,1])
	fig.set_size_inches([3,3])


	pylab.tight_layout()
	if faimsSync:
		path = 'files/app/'
	else:
		path = 'files/server/'

	def update_img(n):    	
		tmp = numpy.random.rand(300,300)
		im.set_data(tmp)
		return im

	#legend(loc=0)
	photoUuid = uuid4()		
	if thumbnail:
		fig.savefig("%s%s_%s.thumbnail.jpg" % (path,photoUuid,filename))
	ani = animation.FuncAnimation(fig,update_img,20,interval=3)
	writer = animation.writers['ffmpeg'](fps=3)

	ani.save('%s%s_%s.mp4' % (path,photoUuid,filename),writer=writer,dpi=dpi)
	pylab.plt.close()
	return "%s%s_%s.mp4" % (path,photoUuid,filename)

def plot_temp():
    data = read_temps()
    dates, values = map(np.array, zip(*[(d['date'], d['temperature'])
                                        for d in data]))
    tmp = (date2num(dates) % 1.0)*24.0
    ii = np.where((tmp > 0) & (tmp < 8))[0]
    continuum = get_continuum(dates, dates[ii], values[ii])
    
    
    setup(figsize=(12,6))
    
    setupplot(subplt=(1,2,1), autoticks=True, xlabel='Date',)
    pylab.plot(dates, values)
    pylab.plot(dates[ii], values[ii], '.r')
    pylab.plot(dates, continuum, '.k')
    plot_weather(np.min(date2num(dates)))
    # pylab.plot(dates, values-continuum+38, '.r')
    dateticks('%Y.%m.%d')
    
    
    setupplot(subplt=(2,2,2), autoticks=False, xlabel='Hour of Day')
    pylab.plot(tmp, values, '.')
    setupplot(subplt=(2,2,2), ylabel='', secondax=True)
    
    setupplot(subplt=(2,2,4), autoticks=False, xlabel='Hour of Day')
    sc = pylab.scatter(tmp, values-continuum+TNORM, 
                       c=date2num(dates)-np.min(date2num(dates)), s=15,
                       marker='.', edgecolor='none',
                       label='Days since Start')
    
    setupplot(subplt=(2,2,4), ylabel='', secondax=True)
    hcolorbar(sc, axes=[0.75, 0.42, 0.1, 0.01])
    
    pylab.tight_layout()
    pylab.show()

def set_font(font):
    pl.xlabel('String length', fontproperties=font)
    pl.tight_layout()
    for label in pl.axes().get_xticklabels():
        label.set_fontproperties(font)
    for label in pl.axes().get_yticklabels():
        label.set_fontproperties(font)

def printHeatMap(marginals, words, outFile):
    N = len(words)
    words_uni = [i.decode('UTF-8') for i in words]
    heatmap = np.zeros((N+1, N+1))
    for chart in marginals:
        heatmap[chart[0], chart[1]] = math.log(marginals[chart])
    fig, ax = plt.subplots()    
    mask = np.tri(heatmap.shape[0], k=0)
    heatmap = np.ma.array(heatmap, mask=mask)
    cmap = plt.cm.get_cmap('RdBu')
    cmap.set_bad('w')
    im = ax.pcolor(heatmap, cmap=cmap, alpha=0.8)
    font = mpl.font_manager.FontProperties(fname='/usr0/home/avneesh/spectral-scfg/data/wqy-microhei.ttf')
    ax.grid(True)
    ax.set_ylim([0,N])
    ax.invert_yaxis()
    ax.set_yticks(np.arange(heatmap.shape[1]-1)+0.5, minor=False)
    ax.set_yticklabels(words_uni, minor=False, fontproperties=font)
    ax.set_xticks(np.arange(heatmap.shape[0])+0.5, minor=True)
    ax.set_xticklabels(np.arange(heatmap.shape[0]), minor=True)
    ax.set_xticks([])
    cbar = fig.colorbar(im, use_gridspec=True)
    cbar.set_label('ln(sum)')
    ax.set_xlabel('Span End')
    ax.xaxis.set_label_position('top')
    ax.xaxis.tick_top()
    plt.ylabel('Span starting at word: ')
    plt.tight_layout()
    #ax.set_title('CKY Heat Map: Node Marginals')
    fig.savefig(outFile)    

def walker_steps(samples, nwalkers, limit):
    ur = N.load('/Users/becky/Projects/Green-Valley-Project/bayesian/find_t_tau_old/colour_plot/ur.npy')
    nuv = N.load('/Users/becky/Projects/Green-Valley-Project/bayesian/find_t_tau_old/colour_plot/nuvu.npy')
    s = samples.reshape(nwalkers, -1, 4)
    s = s[:,:limit,:]
    fig = P.figure(figsize=(9,9))
    ax1 = P.subplot(221, autoscale_on = False, aspect='auto', xlim=[0,13.8], ylim=[0,3])
    ax2 = P.subplot(222,  autoscale_on = False, aspect='auto', xlim=[0,13.8], ylim=[0,3])
    ax3 = P.subplot(223,  autoscale_on = False, aspect='auto', xlim=[0,13.8], ylim=[0,3])
    ax4 = P.subplot(224,  autoscale_on = False, aspect='auto', xlim=[0,13.8], ylim=[0,3])
    ax1.imshow(ur, origin='lower', aspect='auto', extent=[0, 13.8, 0, 3])
    ax2.imshow(ur, origin='lower', aspect='auto', extent=[0, 13.8, 0, 3])
    ax3.imshow(nuv, origin='lower', aspect='auto', extent=[0, 13.8, 0, 3])
    ax4.imshow(nuv, origin='lower', aspect='auto', extent=[0, 13.8, 0, 3])
    for n in range(len(s)):
        ax1.plot(s[n,:,0], s[n,:,1], 'k', alpha=0.5)
        ax2.plot(s[n,:,2], s[n,:,3], 'k', alpha=0.5)
        ax3.plot(s[n,:,0], s[n,:,1], 'k', alpha=0.5)
        ax4.plot(s[n,:,2], s[n,:,3], 'k', alpha=0.5)
    ax1.set_xlabel(r'$t_{bar}$')
    ax1.set_ylabel(r'$\tau_{bar}$')
    ax2.set_xlabel(r'$t_{no bar}$')
    ax2.set_ylabel(r'$\tau_{no bar}$')
    ax3.set_xlabel(r'$t_{bar}$')
    ax3.set_ylabel(r'$\tau_{bar}$')
    ax4.set_xlabel(r'$t_{no bar}$')
    ax4.set_ylabel(r'$\tau_{no bar}$')
    P.tight_layout()
    save_fig = '/Users/becky/Projects/Green-Valley-Project/bayesian/find_t_tau/bars/walkers_2d_steps_bars_'+str(nwalkers)+'_'+str(time.strftime('%H_%M_%d_%m_%y'))+'.pdf'
    fig.savefig(save_fig)
    return fig

def pca(inF,MIN):
    df = pd.read_table(inF, header=0)
    dc = list(df.columns)
    dc[0]='GeneID'
    df.columns = dc
    
    print(df.shape)
    sel = ~((df.ix[:,2] < MIN) & (df.ix[:,3]< MIN) & (df.ix[:,4]< MIN) & (df.ix[:,5]< MIN) & (df.ix[:,6]< MIN) & (df.ix[:,7]< MIN) & (df.ix[:,8]< MIN) & (df.ix[:,9]< MIN))
    df = df.ix[sel,:]
    print(df.shape)
    
    X = df.ix[:,2:df.shape[1]].values.T
    y = df.columns[2:df.shape[1]].values
    X_std = StandardScaler().fit_transform(X)
    
    #pca = PCA(n_components=2)
    pca = PCA()
    Y_sklearn = pca.fit_transform(X_std)
    
    
    fig = plt.figure()
    plt.style.use('ggplot')
    #plt.style.use('seaborn-whitegrid')
    ax = fig.add_subplot(111)
    for lab, col in zip(y,('red','red', 'green','green', 'blue','blue','m','m')):
        ax.scatter(Y_sklearn[y==lab, 0],Y_sklearn[y==lab, 1],label=lab,c=col, s=80)
    
    
    ax.set_xlabel('Principal Component 1 : %.2f'%(pca.explained_variance_ratio_[0]*100) + '%')
    ax.set_ylabel('Principal Component 2 : %.2f'%(pca.explained_variance_ratio_[1]*100) + '%')
    ax.legend(loc='lower right', prop={'size':8})
    plt.tight_layout()
    plt.savefig(inF + '-RNASeq-MIN' + str(MIN) + '.pdf')

def corner_plot(s, labels):
    x, y = s[:,0], s[:,1]
    fig = P.figure(figsize=(10,10))
    ax2 = P.subplot(223)
    ax2.set_xlabel(labels[0])
    ax2.set_ylabel(labels[1])
    im = triangle.histo2d(x, y, ax=ax2, extent=[[0, 13.807108309208775],[0, 3.0]])
    [l.set_rotation(45) for l in ax2.get_xticklabels()]
    [j.set_rotation(45) for j in ax2.get_yticklabels()]
    ax1 = P.subplot(221, xlim=[0, 13.807108309208775])
    ax1.tick_params(axis='x', labelbottom='off')
    ax1.tick_params(axis='y', labelleft='off')
    ax1.hist(x, bins=50, histtype='step', color='k', range=(0, 13.807108309208775))
    ax3 = P.subplot(224)
    ax3.tick_params(axis='x', labelbottom='off')
    ax3.tick_params(axis='y', labelleft='off')
    ax3.hist(y, bins=50, orientation='horizontal', histtype='step',color='k', range=(0,3))
    P.subplots_adjust(wspace=0.05)
    P.subplots_adjust(hspace=0.05)
    cbar_ax = fig.add_axes([0.55, 0.565, 0.02, 0.405])
    cb = fig.colorbar(im, cax = cbar_ax)
    cb.solids.set_edgecolor('face')
    cb.set_label(r'predicted SFR $[M_{\odot} yr^{-1}]$', labelpad = 20, fontsize=16)
    P.tight_layout()
    return fig

def plot_file_color(base, thin=True, start=0, size=14, save=False):
    conf, track, pegs = load(base)

    fig = pl.figure(figsize=(size,size*conf['top']/conf['wall']))

    track = track[start:]
    x = track[:,0];   y = track[:,1]
    t = np.linspace(0,1,x.shape[0])
    points = np.array([x,y]).transpose().reshape(-1,1,2)
    segs = np.concatenate([points[:-1],points[1:]],axis=1)
    lc = LineCollection(segs, linewidths=0.25, cmap=pl.cm.coolwarm)
    lc.set_array(t)
    pl.gca().add_collection(lc)

    #pl.scatter(x, y, c=np.arange(len(x)),linestyle='-',cmap=pl.cm.coolwarm)
    #pl.plot(track[-1000000:,0], track[-1000000:,1], '-', linewidth=0.0125, alpha=0.8)
    for peg in pegs:
        pl.gca().add_artist(pl.Circle(peg, conf['radius'], color='k', alpha=0.3))
    pl.xlim(0, conf['wall'])
    pl.ylim(0, conf['top'])
    pl.xticks([])
    pl.yticks([])
    pl.tight_layout()
    pl.show()
    if save:
        pl.savefig(base+".png", dpi=200)

    def plot(self):
        f = pylab.figure(figsize=(8,4))
        co = [] #colors container
        for zScore, r in itertools.izip(self.zScores, self.log2Ratio):
            if zScore < self.pCut:
                if r > 0:
                    co.append(Colors().greenColor)
                elif r < 0:
                    co.append(Colors().redColor)
                else:
                    raise Exception
            else:
                co.append(Colors().blueColor)

        #print "Probability this is from a normal distribution: %.3e" %stats.normaltest(self.log2Ratio)[1]
        ax = f.add_subplot(121)
        pylab.axvline(self.meanLog2Ratio, color=Colors().redColor)
        pylab.axvspan(self.meanLog2Ratio-(2*self.stdLog2Ratio), 
                      self.meanLog2Ratio+(2*self.stdLog2Ratio), color=Colors().blueColor, alpha=0.2)
        his = pylab.hist(self.log2Ratio, bins=50, color=Colors().blueColor)
        pylab.xlabel("log2 Ratio %s/%s" %(self.sampleNames[1], self.sampleNames[0]))
        pylab.ylabel("Frequency")
        
        ax = f.add_subplot(122, aspect='equal')
        pylab.scatter(self.genes1, self.genes2, c=co, alpha=0.5)        
        pylab.ylabel("%s RPKM" %self.sampleNames[1])
        pylab.xlabel("%s RPKM" %self.sampleNames[0])
        pylab.yscale('log')
        pylab.xscale('log')
        pylab.tight_layout()