Python pylab.colorbar() Examples

def plot_confusion_matrix(self, matrix, labels):
        if not self.to_save and not self.to_show:
            return

        pylab.figure()
        pylab.imshow(matrix, interpolation='nearest', cmap=pylab.cm.jet)
        pylab.title("Confusion Matrix")

        for i, vi in enumerate(matrix):
            for j, vj in enumerate(vi):
                pylab.annotate("%.1f" % vj, xy=(j, i), horizontalalignment='center', verticalalignment='center', fontsize=9)

        pylab.colorbar()

        classes = np.arange(len(labels))
        pylab.xticks(classes, labels)
        pylab.yticks(classes, labels)

        pylab.ylabel('Expected label')
        pylab.xlabel('Predicted label') 

def plot_wt_layout(wt_layout, borders=None, depth=None):
    fig = plt.figure(figsize=(6,6), dpi=2000)
    fs = 14
    ax = plt.subplot(111)

    if depth is not None:
        N = 100
        X, Y = plt.meshgrid(plt.linspace(depth[:,0].min(), depth[:,0].max(), N), 
                            plt.linspace(depth[:,1].min(), depth[:,1].max(), N))
        Z = plt.griddata(depth[:,0],depth[:,1],depth[:,2],X,Y, interp='linear')
        plt.contourf(X,Y,Z, label='depth [m]')
        plt.colorbar().set_label('water depth [m]')
    #ax.plot(wt_layout.wt_positions[:,0], wt_layout.wt_positions[:,1], 'or', label='baseline position')
    
    ax.scatter(wt_layout.wt_positions[:,0], wt_layout.wt_positions[:,1], wt_layout._wt_list('rotor_diameter'), label='baseline position')

    if borders is not None:
        ax.plot(borders[:,0], borders[:,1], 'r--', label='border')
        
    ax.set_xlabel('x [m]'); 
    ax.set_ylabel('y [m]')
    ax.axis('equal');
    ax.legend(loc='lower left') 

def plot_functional_map(C, newfig=True):
    vmax = max(np.abs(C.max()), np.abs(C.min()))
    vmin = -vmax
    C = ((C - vmin) / (vmax - vmin)) * 2 - 1
    if newfig:
        pl.figure(figsize=(5,5))
    else:
        pl.clf()
    ax = pl.gca()
    pl.pcolor(C[::-1], edgecolor=(0.9, 0.9, 0.9, 1), lw=0.5,
              vmin=-1, vmax=1, cmap=nice_mpl_color_map())
    # colorbar
    tick_locs   = [-1., 0.0, 1.0]
    tick_labels = ['min', 0, 'max']
    bar = pl.colorbar()
    bar.locator = matplotlib.ticker.FixedLocator(tick_locs)
    bar.formatter = matplotlib.ticker.FixedFormatter(tick_labels)
    bar.update_ticks()
    ax.set_aspect(1)
    pl.xticks([])
    pl.yticks([])
    if newfig:
        pl.show() 

def heatmap(df,fname=None,cmap='seismic',log=False):
    """Plot a heat map"""

    from matplotlib.colors import LogNorm
    f=plt.figure(figsize=(8,8))
    ax=f.add_subplot(111)
    norm=None
    df=df.replace(0,.1)
    if log==True:
        norm=LogNorm(vmin=df.min().min(), vmax=df.max().max())
    hm = ax.pcolor(df,cmap=cmap,norm=norm)
    plt.colorbar(hm,ax=ax,shrink=0.6,norm=norm)
    plt.yticks(np.arange(0.5, len(df.index), 1), df.index)
    plt.xticks(np.arange(0.5, len(df.columns), 1), df.columns, rotation=90)
    #ax.axvline(4, color='gray'); ax.axvline(8, color='gray')
    plt.tight_layout()
    if fname != None:
        f.savefig(fname+'.png')
    return ax 

def plotOutput(layer,feed_dict,fieldShape=None,channel=None,figOffset=1,cmap=None):
	# Output summary
	try:
		W = layer.output
	except:
		W = layer
	wp = W.eval(feed_dict=feed_dict);
	if len(np.shape(wp)) < 4:		# Fully connected layer, has no shape
		temp = np.zeros(np.product(fieldShape)); temp[0:np.shape(wp.ravel())[0]] = wp.ravel()
		fields = np.reshape(temp,[1]+fieldShape)
	else:			# Convolutional layer already has shape
		wp = np.rollaxis(wp,3,0)
		features, channels, iy,ix = np.shape(wp)
		if channel is not None:
			fields = wp[:,channel,:,:]
		else:
			fields = np.reshape(wp,[features*channels,iy,ix])

	perRow = int(math.floor(math.sqrt(fields.shape[0])))
	perColumn = int(math.ceil(fields.shape[0]/float(perRow)))
	fields2 = np.vstack([fields,np.zeros([perRow*perColumn-fields.shape[0]] + list(fields.shape[1:]))])
	tiled = []
	for i in range(0,perColumn*perRow,perColumn):
		tiled.append(np.hstack(fields2[i:i+perColumn]))

	tiled = np.vstack(tiled)
	if figOffset is not None:
		mpl.figure(figOffset); mpl.clf(); 

	mpl.imshow(tiled,cmap=cmap); mpl.title('%s Output' % layer.name); mpl.colorbar(); 

def plot_confusion_matrix(cm, title='Confusion matrix', cmap=plt.cm.Blues):
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(iris.target_names))
    plt.xticks(tick_marks, iris.target_names, rotation=45)
    plt.yticks(tick_marks, iris.target_names)
    plt.tight_layout()
    plt.ylabel('True label')
    plt.xlabel('Predicted label') 

def nice_imshow(ax, data, vmin=None, vmax=None, cmap=None):
    """Wrapper around pl.imshow"""
    if cmap is None:
        cmap = cm.jet
    if vmin is None:
        vmin = data.min()
    if vmax is None:
        vmax = data.max()
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    im = ax.imshow(data, vmin=vmin, vmax=vmax, interpolation='nearest', cmap=cmap)
    pl.colorbar(im, cax=cax) 

def plotHeatMap(data, col='KL without null', label=''):

    #Compute and collate medians
    sel_cols = [x for x in data.columns if col in x]
    cmp_meds = data[sel_cols].median(axis=0)
    samples = sortSampleNames(getUniqueSamples(sel_cols))
    cell_lines = ['CHO', 'E14TG2A', 'BOB','RPE1', 'HAP1','K562','eCAS9','TREX2']
    sample_idxs = [(cell_lines.index(parseSampleName(x)[0]),x) for x in getUniqueSamples(sel_cols)]
    sample_idxs.sort()
    samples = [x[1] for x in sample_idxs]

    N = len(samples)
    meds = np.zeros((N,N))
    for colname in sel_cols:
        dir1, dir2 = getDirsFromFilename(colname.split('$')[-1])
        idx1, idx2 = samples.index(dir1), samples.index(dir2)
        meds[idx1,idx2] = cmp_meds[colname]
        meds[idx2,idx1] = cmp_meds[colname]

    for i in range(N):
        print(' '.join(['%.2f' % x for x in meds[i,:]]))
    print( np.median(meds[:,:-4],axis=0))

	#Display in Heatmap
    PL.figure(figsize=(5,5))
    PL.imshow(meds, cmap='hot_r', vmin = 0.0, vmax = 3.0, interpolation='nearest')
    PL.colorbar()
    PL.xticks(range(N))
    PL.yticks(range(N))
    PL.title("Median KL") # between %d mutational profiles (for %s with >%d mutated reads)" % (col, len(data), label, MIN_READS))
    ax1 = PL.gca()
    ax1.set_yticklabels([getSimpleName(x) for x in samples], rotation='horizontal')
    ax1.set_xticklabels([getSimpleName(x) for x in samples], rotation='vertical')
    PL.subplots_adjust(left=0.25,right=0.95,top=0.95, bottom=0.25)
    PL.show(block=False) 
    saveFig('median_kl_heatmap_cell_lines') 

def dendogram(correlation_matrix,
	save_as = '',
	figsize=(50,50),
	):

	correlation_matrix = path.abspath(path.expanduser(correlation_matrix))

	y = hier_clustering(correlation_matrix, method='centroid')
	z = hier_clustering(y, orientation='right')

	fig = pylab.figure(figsize=figsize)
	ax_1 = fig.add_axes([0.1,0.1,0.2,0.8])
	ax_1.set_xticks([])
	ax_1.set_yticks([])

	ax_2 = fig.add_axes([0.3,0.1,0.6,0.8])
	index = z['leaves']
	correlation_matrix = correlation_matrix[index,:]
	correlation_matrix = correlation_matrix[:,index]
	im = ax_2.matshow(correlation_matrix, aspect='auto', origin='lower')
	ax_2.set_xticks([])
	ax_2.set_yticks([])

	ax_color = fig.add_axes([0.91,0.1,0.02,0.8])
	colorbar = pylab.colorbar(im, cax=ax_color)
	colorbar.ax.tick_params(labelsize=75)

	# Display and save figure.
	if(save_as):
		fig.savefig(path.abspath(path.expanduser(save_as))) 

def examplify(cls, fname):
        """
        example of how to use
        """
        logger.info("examplify starts")
        
        # model
        model = MultiClassLogisticRegression(epsilon=0.01, n_scan = 100)
        
        # learn
        st = time.time()
        model.learn(fname)
        et = time.time()
        print "learning time: %d [s]" % ((et - st)/1000)

        # predict
        y_pred = model.predict(fname)
        
        # confusion matrix
        y_label = np.loadtxt(fname, delimiter=" ")[:, 0]

        cm = confusion_matrix(y_label, y_pred)
        #pl.matshow(cm)
        #pl.title('Confusion matrix')
        #pl.colorbar()
        #pl.ylabel('True label')
        #pl.xlabel('Predicted label')
        #pl.show()

        print cm
        print "accurary: %d [%%]" % (np.sum(cm.diagonal()) * 100.0/np.sum(cm))
        logger.info("examplify finished") 

def plot(self, filename=None, vmin=None, vmax=None, cmap='jet_r'):
        import pylab
        pylab.clf()
        pylab.imshow(-np.log10(self.results[self._start_y:,:]), 
            origin="lower",
            aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax)
        pylab.colorbar()

        # Fix xticks
        XMAX = float(self.results.shape[1])  # The max integer on xaxis
        xpos = list(range(0, int(XMAX), int(XMAX/5)))
        xx = [int(this*100)/100 for this in np.array(xpos) / XMAX * self.duration]
        pylab.xticks(xpos, xx, fontsize=16)

        # Fix yticks
        YMAX = float(self.results.shape[0])  # The max integer on xaxis
        ypos = list(range(0, int(YMAX), int(YMAX/5)))
        yy = [int(this) for this in np.array(ypos) / YMAX * self.sampling]
        pylab.yticks(ypos, yy, fontsize=16)

        #pylab.yticks([1000,2000,3000,4000], [5500,11000,16500,22000], fontsize=16)
        #pylab.title("%s echoes" %  filename.replace(".png", ""), fontsize=25)
        pylab.xlabel("Time (seconds)", fontsize=25)
        pylab.ylabel("Frequence (Hz)", fontsize=25)
        pylab.tight_layout()
        if filename:
            pylab.savefig(filename) 

def contour_plot(func):
    rose = func()
    XS, YS = plt.meshgrid(np.linspace(-2, 2, 20), np.linspace(-2,2, 20));
    ZS = np.array([rose(x1=x, x2=y).f_xy for x,y in zip(XS.flatten(),YS.flatten())]).reshape(XS.shape);
    plt.contourf(XS, YS, ZS, 50);
    plt.colorbar() 

def plotOutput(layer,feed_dict,fieldShape=None,channel=None,figOffset=1,cmap=None):
	# Output summary
	W = layer.output
	wp = W.eval(feed_dict=feed_dict);
	if len(np.shape(wp)) < 4:		# Fully connected layer, has no shape
		temp = np.zeros(np.product(fieldShape)); temp[0:np.shape(wp.ravel())[0]] = wp.ravel()
		fields = np.reshape(temp,[1]+fieldShape)
	else:			# Convolutional layer already has shape
		wp = np.rollaxis(wp,3,0)
		features, channels, iy,ix = np.shape(wp)
		if channel is not None:
			fields = wp[:,channel,:,:]
		else:
			fields = np.reshape(wp,[features*channels,iy,ix])

	perRow = int(math.floor(math.sqrt(fields.shape[0])))
	perColumn = int(math.ceil(fields.shape[0]/float(perRow)))
	fields2 = np.vstack([fields,np.zeros([perRow*perColumn-fields.shape[0]] + list(fields.shape[1:]))])
	tiled = []
	for i in range(0,perColumn*perRow,perColumn):
		tiled.append(np.hstack(fields2[i:i+perColumn]))

	tiled = np.vstack(tiled)
	if figOffset is not None:
		mpl.figure(figOffset); mpl.clf();

	mpl.imshow(tiled,cmap=cmap); mpl.title('%s Output' % layer.name); mpl.colorbar(); 

def plotFields(layer,fieldShape=None,channel=None,maxFields=25,figName='ReceptiveFields',cmap=None,padding=0.01):
	# Receptive Fields Summary
	W = layer.W
	wp = W.eval().transpose();
	if len(np.shape(wp)) < 4:		# Fully connected layer, has no shape
		fields = np.reshape(wp,list(wp.shape[0:-1])+fieldShape)
	else:			# Convolutional layer already has shape
		features, channels, iy, ix = np.shape(wp)
		if channel is not None:
			fields = wp[:,channel,:,:]
		else:
			fields = np.reshape(wp,[features*channels,iy,ix])

	fieldsN = min(fields.shape[0],maxFields)
	perRow = int(math.floor(math.sqrt(fieldsN)))
	perColumn = int(math.ceil(fieldsN/float(perRow)))

	fig = mpl.figure(figName); mpl.clf()

	# Using image grid
	from mpl_toolkits.axes_grid1 import ImageGrid
	grid = ImageGrid(fig,111,nrows_ncols=(perRow,perColumn),axes_pad=padding,cbar_mode='single')
	for i in range(0,fieldsN):
		im = grid[i].imshow(fields[i],cmap=cmap);

	grid.cbar_axes[0].colorbar(im)
	mpl.title('%s Receptive Fields' % layer.name)

	# old way
	# fields2 = np.vstack([fields,np.zeros([perRow*perColumn-fields.shape[0]] + list(fields.shape[1:]))])
	# tiled = []
	# for i in range(0,perColumn*perRow,perColumn):
	# 	tiled.append(np.hstack(fields2[i:i+perColumn]))
	#
	# tiled = np.vstack(tiled)
	# mpl.figure(figOffset); mpl.clf(); mpl.imshow(tiled,cmap=cmap); mpl.title('%s Receptive Fields' % layer.name); mpl.colorbar();
	mpl.figure(figName+' Total'); mpl.clf(); mpl.imshow(np.sum(np.abs(fields),0),cmap=cmap); mpl.title('%s Total Absolute Input Dependency' % layer.name); mpl.colorbar() 

def plotFields(layer,fieldShape=None,channel=None,figOffset=1,cmap=None,padding=0.01):
	# Receptive Fields Summary
	try:
		W = layer.W
	except:
		W = layer
	wp = W.eval().transpose();
	if len(np.shape(wp)) < 4:		# Fully connected layer, has no shape
		fields = np.reshape(wp,list(wp.shape[0:-1])+fieldShape)	
	else:			# Convolutional layer already has shape
		features, channels, iy, ix = np.shape(wp)
		if channel is not None:
			fields = wp[:,channel,:,:]
		else:
			fields = np.reshape(wp,[features*channels,iy,ix])

	perRow = int(math.floor(math.sqrt(fields.shape[0])))
	perColumn = int(math.ceil(fields.shape[0]/float(perRow)))

	fig = mpl.figure(figOffset); mpl.clf()
	
	# Using image grid
	from mpl_toolkits.axes_grid1 import ImageGrid
	grid = ImageGrid(fig,111,nrows_ncols=(perRow,perColumn),axes_pad=padding,cbar_mode='single')
	for i in range(0,np.shape(fields)[0]):
		im = grid[i].imshow(fields[i],cmap=cmap); 

	grid.cbar_axes[0].colorbar(im)
	mpl.title('%s Receptive Fields' % layer.name)
	
	# old way
	# fields2 = np.vstack([fields,np.zeros([perRow*perColumn-fields.shape[0]] + list(fields.shape[1:]))])
	# tiled = []
	# for i in range(0,perColumn*perRow,perColumn):
	# 	tiled.append(np.hstack(fields2[i:i+perColumn]))
	# 
	# tiled = np.vstack(tiled)
	# mpl.figure(figOffset); mpl.clf(); mpl.imshow(tiled,cmap=cmap); mpl.title('%s Receptive Fields' % layer.name); mpl.colorbar();
	mpl.figure(figOffset+1); mpl.clf(); mpl.imshow(np.sum(np.abs(fields),0),cmap=cmap); mpl.title('%s Total Absolute Input Dependency' % layer.name); mpl.colorbar() 

def hierarchical_clustering(mat, method='average', cluster_distance=True, labels=None, thres=0.65):
    """
    Performs hierarchical clustering based on distance matrix 'mat' using the method specified by 'method'.
    Optional argument 'labels' may specify a list of labels. If cluster_distance is True, the clustering is
    performed on the distance matrix using euclidean distance. Otherwise, mat specifies the distance matrix for
    clustering. Adapted from
    http://stackoverflow.com/questions/7664826/how-to-get-flat-clustering-corresponding-to-color-clusters-in-the-dendrogram-cre
    Not subjected to copyright.
    """
    D = numpy.array(mat)
    if not cluster_distance:
        Dtriangle = scipy.spatial.distance.squareform(D)
    else:
        Dtriangle = scipy.spatial.distance.pdist(D, metric='euclidean')
    fig = pylab.figure(figsize=(8, 8))
    ax1 = fig.add_axes([0.09, 0.1, 0.2, 0.6])
    Y = sch.linkage(Dtriangle, method=method)
    Z1 = sch.dendrogram(Y, orientation='right', color_threshold=thres*max(Y[:, 2]))
    ax1.set_xticks([])
    ax1.set_yticks([])
    ax2 = fig.add_axes([0.3, 0.71, 0.6, 0.2])
    Y = sch.linkage(Dtriangle, method=method)
    Z2 = sch.dendrogram(Y, color_threshold=thres*max(Y[:, 2]))
    ax2.set_xticks([])
    ax2.set_yticks([])
    axmatrix = fig.add_axes([0.3, 0.1, 0.6, 0.6])
    idx1 = Z1['leaves']
    idx2 = Z2['leaves']
    D = D[idx1, :]
    D = D[:, idx2]
    im = axmatrix.matshow(D, aspect='auto', origin='lower', cmap=pylab.get_cmap('jet_r'))
    if labels is None:
        axmatrix.set_xticks([])
        axmatrix.set_yticks([])
    else:
        axmatrix.set_xticks(range(len(labels)))
        lab = [labels[idx1[m]] for m in range(len(labels))]
        axmatrix.set_xticklabels(lab)
        axmatrix.set_yticks(range(len(labels)))
        axmatrix.set_yticklabels(lab)
        for tick in pylab.gca().xaxis.iter_ticks():
            tick[0].label2On = False
            tick[0].label1On = True
            tick[0].label1.set_rotation('vertical')
        for tick in pylab.gca().yaxis.iter_ticks():
            tick[0].label2On = True
            tick[0].label1On = False
    axcolor = fig.add_axes([0.91, 0.1, 0.02, 0.6])
    pylab.colorbar(im, cax=axcolor)
    pylab.show()
    return Z1 

def result_analysis(prediction, train_truth, valid_truth, test_truth, 
                    verbose=False):
    assert prediction.shape == test_truth.shape
    print "Detailed information in each category:"
    print "                          Number of Samples"
    print "Class No.       TRAIN   VALID   TEST    RightCount   RightRate"
    for i in xrange(test_truth.min(), test_truth.max()+1):
        right_prediction = ( (test_truth-prediction) == 0 )
        right_count = numpy.sum(((test_truth==i) * right_prediction)*1)
        print "%d\t\t%d\t%d\t%d\t%d\t%f" % \
            (i, 
             numpy.sum((train_truth==i)*1), 
             numpy.sum((valid_truth==i)*1), 
             numpy.sum((test_truth==i)*1),
             right_count,
             right_count * 1.0 / numpy.sum((test_truth==i)*1)
            )
    
    total_right_count = numpy.sum(right_prediction*1)
    print "Overall\t\t%d\t%d\t%d\t%d\t%f" % \
            (train_truth.size, 
             valid_truth.size, 
             test_truth.size, 
             total_right_count,
             total_right_count * 1.0 / test_truth.size
            )
    
    cm = confusion_matrix(test_truth, prediction)
    pr_a = cm.trace()*1.0 / test_truth.size
    pr_e = ((cm.sum(axis=0)*1.0/test_truth.size) * \
            (cm.sum(axis=1)*1.0/test_truth.size)).sum()
    k = (pr_a - pr_e) / (1 - pr_e)
    print "kappa index of agreement: %f" % k
    print "confusion matrix:"
    print cm
    # Show confusion matrix
    pl.matshow(cm)
    pl.title('Confusion matrix')
    pl.colorbar()
    if verbose:
        pl.show()
    else:
        filename = 'conf_mtx_' + str(time.time()) + '.png'
        pl.savefig(filename)


#------------------------------------------------------------------------------- 

def callback_gppvae0(epoch, history, covs, imgs, ffile):

    # init fig
    pl.figure(1, figsize=(8, 8))
    pl.subplot(4, 4, 1)
    pl.title("loss")
    pl.plot(history["loss"])
    pl.subplot(4, 4, 2)
    pl.title("vars")
    pl.plot(sp.array(history["vs"])[:, 0], "r")
    pl.plot(sp.array(history["vs"])[:, 1], "k")
    pl.subplot(4, 4, 5)
    pl.title("mse_out")
    pl.plot(history["mse_out"])

    pl.subplot(4, 4, 9)
    pl.title("XX")
    pl.imshow(covs["XX"], vmin=-0.4, vmax=1)
    pl.colorbar()
    pl.subplot(4, 4, 10)
    pl.title("WW")
    pl.imshow(covs["WW"], vmin=-0.4, vmax=1)
    pl.colorbar()

    Yv, Rv = imgs["Yv"], imgs["Yo"]

    # make plot
    pl.subplot(4, 2, 2)
    _img = _compose(Yv[0:6], Rv[0:6])
    pl.imshow(_img)

    pl.subplot(4, 2, 4)
    _img = _compose(Yv[6:12], Rv[6:12])
    pl.imshow(_img)

    pl.subplot(4, 2, 6)
    _img = _compose(Yv[12:18], Rv[12:18])
    pl.imshow(_img)

    pl.subplot(4, 2, 8)
    _img = _compose(Yv[18:24], Rv[18:24])
    pl.imshow(_img)

    pl.savefig(ffile)
    pl.close() 

def callback_gppvae(epoch, history, covs, imgs, ffile):

    # init fig
    pl.figure(1, figsize=(8, 8))
    pl.subplot(4, 4, 1)
    pl.title("loss")
    pl.plot(history["loss"], "k")
    pl.subplot(4, 4, 2)
    pl.title("vars")
    pl.plot(sp.array(history["vs"])[:, 0], "r")
    pl.plot(sp.array(history["vs"])[:, 1], "k")
    pl.ylim(0, 1)
    pl.subplot(4, 4, 5)
    pl.title("recon_term")
    pl.plot(history["recon_term"], "k")
    pl.subplot(4, 4, 6)
    pl.title("gp_nll")
    pl.plot(history["gp_nll"], "k")
    pl.subplot(4, 4, 9)
    pl.title("mse_out")
    pl.plot(history["mse_out"], "k")
    pl.ylim(0, 0.1)
    pl.subplot(4, 4, 10)
    pl.title("mse")
    pl.plot(history["mse"], "k")
    pl.plot(history["mse_val"], "r")
    pl.ylim(0, 0.01)

    pl.subplot(4, 4, 13)
    pl.title("XX")
    pl.imshow(covs["XX"], vmin=-0.4, vmax=1)
    pl.colorbar()
    pl.subplot(4, 4, 14)
    pl.title("WW")
    pl.imshow(covs["WW"], vmin=-0.4, vmax=1)
    pl.colorbar()

    Yv, Yr, Rv = imgs["Yv"], imgs["Yr"], imgs["Yo"]

    # make plot
    pl.subplot(4, 2, 2)
    _img = _compose_multi([Yv[0:6], Yr[0:6], Rv[0:6]])
    pl.imshow(_img)

    pl.subplot(4, 2, 4)
    _img = _compose_multi([Yv[6:12], Yr[6:12], Rv[6:12]])
    pl.imshow(_img)

    pl.subplot(4, 2, 6)
    _img = _compose_multi([Yv[12:18], Yr[12:18], Rv[12:18]])
    pl.imshow(_img)

    pl.subplot(4, 2, 8)
    _img = _compose_multi([Yv[18:24], Yr[18:24], Rv[18:24]])
    pl.imshow(_img)

    pl.savefig(ffile)
    pl.close() 

def gs_mod(idata,itera=10,osize=256):
    """Modiffied Gerchberg-Saxton algorithm to calculate DOEs
    
    Calculates the phase distribution in a object plane to obtain an 
    specific amplitude distribution in the target plane. It uses a 
    FFT to calculate the field propagation.
    The wavefront at the DOE plane is assumed as a plane wave.
    This algorithm leaves a window around the image plane to allow the 
    noise to move there. It only optimises the center of the image.
    
    **ARGUMENTS:**
    
        ========== ======================================================
        idata      numpy array containing the target amplitude distribution 
        itera      Maximum number of iterations
        osize      Size of the center of the image to be optimized
                   It should be smaller than the image itself.
        ========== ======================================================
    """
    M,N=idata.shape
    cut=osize//2
    
    
    zone=zeros_like(idata)
    zone[M/2-cut:M/2+cut,N/2-cut:N/2+cut]=1
    zone=zone.astype(bool)

    mask=exp(2.j*pi*random(idata.shape))
    mask[zone]=0
    
    #~ imshow(abs(mask)),colorbar()
    
    fdata=fftshift(fft2(ifftshift(idata+mask))) #Nota, colocar esta mascara es muy importante, por que si no  no converge tan rapido
    
    e=1000
    ea=1000
    for i in range (itera):
        fdata=exp(1.j*angle(fdata))

        rdata=ifftshift(ifft2(fftshift(fdata)))
        #~ e= (abs(rdata[zone])-idata[zone]).std()
        #~ if e>ea: 
           #~ 
            #~ break
        ea=e
        rdata[zone]=exp(1.j*angle(rdata[zone]))*(idata[zone])        
        fdata=fftshift(fft2(ifftshift(rdata)))   
    fdata=exp(1.j*angle(fdata))
    return fdata 

def plotCorrelations(all_result_outputs, label='', data_label='', y_label='', plot_label='', plot_scatters=False, oligo_id_str='Oligo ID', val_str = 'Cut Rate', total_reads_str= 'Total Reads', scatter_samples={}, sdims=(0,0), scatter_fig=None, add_leg=True):
    
    datas = [x[0][data_label][0] for x in all_result_outputs]
    sample_names = [shortDirLabel(x[1]) for x in all_result_outputs]
    
    merged_data = pd.merge(datas[0],datas[1],how='inner',on=oligo_id_str, suffixes=['', ' 2'])
    for i, data in enumerate(datas[2:]):
        merged_data = pd.merge(merged_data, data,how='inner',on=oligo_id_str, suffixes=['', ' %d' % (i+3)])
    suffix = lambda i: ' %d' % (i+1) if i > 0 else ''

    N = len(sample_names)
    if plot_scatters: 
        if scatter_fig is None: PL.figure()
        else: PL.figure(scatter_fig.number)
        s_dims = (N,N) if len(scatter_samples) == 0 else sdims
    pcorrs, scorrs = np.zeros((N,N)), np.zeros((N,N))
    for i,label1 in enumerate(sample_names):
        for j,label2 in enumerate(sample_names):
            dvs1, dvs2, ids = merged_data[val_str + suffix(i)], merged_data[val_str + suffix(j)], merged_data[oligo_id_str]
            pcorrs[i,j] = pearsonr(dvs1, dvs2)[0]
            scorrs[i,j] = spearmanr(dvs1, dvs2)[0]
            if plot_scatters:
                if (label1, label2) in scatter_samples: idx = scatter_samples[(label1, label2)]
                elif len(scatter_samples) == 0: idx = i*N+j+1
                else: continue
                PL.subplot(s_dims[0],s_dims[1],idx)
                trs1, trs2 = merged_data[total_reads_str + suffix(i)], merged_data[total_reads_str + suffix(j)]
                for thr in [20,50,100,500,1000]:
                    thr_dvs = [(dv1, dv2,id) for (dv1, dv2, tr1, tr2,id) in zip(dvs1,dvs2,trs1,trs2,ids) if (tr1 >= thr and tr2 >= thr) ]
                    pcorrs[i,j] = pearsonr([x[0] for x in thr_dvs], [x[1] for x in thr_dvs])[0]
                    PL.plot([x[0] for x in thr_dvs],[x[1] for x in thr_dvs],'.', label='>%d Reads' % thr)

                PL.plot([0,100],[0,100],'k--')
                PL.xlabel('K562 Replicate A')
                PL.ylabel('K562 Replicate B')
                if add_leg: PL.legend()
                PL.title('%s (%.2f)' % (y_label, pcorrs[i,j]))
    if plot_scatters: 
        PL.subplots_adjust(left=0.05,right=0.95,top=0.9, bottom=0.1, hspace=0.4)
        PL.show(block=False)
        saveFig(plot_label)

    if not plot_scatters:
        PL.figure(figsize=(10,10))
        #PL.subplot(1,2,1)
        PL.imshow(pcorrs, cmap='hot', vmin = 0.0, vmax = 1.0, interpolation='nearest')
        PL.xticks(range(N),sample_names, rotation='vertical')
        PL.yticks(range(N),sample_names)
        PL.title(y_label + ': Pearson')
        PL.colorbar()
        #PL.subplot(1,2,2)
        #PL.imshow(scorrs, cmap='hot', vmin = 0.0, vmax = 1.0, interpolation='nearest')
        #PL.xticks(range(N),sample_names, rotation='vertical')
        #PL.yticks(range(N),sample_names)
        #PL.title(y_label + ': Spearman')
        #PL.colorbar()
        PL.show(block=False)
        PL.savefig(getPlotDir() + '/%s_%s.png' % (plot_label, sanitizeLabel(label)), bbox_inches='tight') 

def plot_glm_results(axis, results, best_rmse, cb):
    axis.cla()
    axis.set_xscale('log')
    axis.set_xlim([1e2, 1e9])
    axis.set_ylim([-0.12, 1.12])
    axis.set_yticks([x / 7. for x in range(0, 8)])
    axis.set_ylabel('Parameter 1:  ' + r'$\alpha$', fontsize=16)
    axis.set_xlabel('Parameter 2:  ' + r'$\lambda$', fontsize=16)
    num_models = min(4000, int(4000 * results.shape[0] / 2570))
    axis.set_title(
        'Elastic Net Models Trained and Evaluated: ' + str(num_models),
        fontsize=16)

    try:
        import seaborn as sns
        sns.set_style("whitegrid")
        import pylab as pl
        from matplotlib.colors import ListedColormap
        cm = ListedColormap(sns.color_palette("RdYlGn", 10).as_hex())
        cf = axis.scatter(
            results['lambda'],
            results['alpha_prime'],
            c=results['rel_acc'],
            cmap=cm,
            vmin=0,
            vmax=1,
            s=60,
            lw=0)
        axis.plot(
            best_rmse['lambda'],
            best_rmse['alpha_prime'],
            'o',
            ms=15,
            mec='k',
            mfc='none',
            mew=2)

        if not cb:
            cb = pl.colorbar(cf, ax=axis)
            cb.set_label(
                'Relative  Validation  Accuracy',
                rotation=270,
                labelpad=18,
                fontsize=16)
        cb.update_normal(cf)
    except:
        #print("plot_glm_results exception -- no frame")
        pass 

def visualize_feature_maps(
    X, anno=[], keypoints=[], axes=[2, 0], save_filename=None, divide=False
):
    nc = np.ceil(np.sqrt(X.shape[2]))  # column
    nr = np.ceil(X.shape[2] / nc)  # row
    nc = int(nc)
    nr = int(nr)
    plt.figure(figsize=(64, 64))
    for i in range(X.shape[2]):
        ax = plt.subplot(nr, nc, i + 1)
        ax.imshow(X[:, :, i], cmap="jet")
        # plt.colorbar()
        for obj in anno:
            draw_box(ax, obj, axes=axes, color="r")

            # plot head orientation
            center_bottom_forward = np.mean(obj.T[[0, 1, 5, 4]], axis=0)
            center_bottom = np.mean(obj.T[[1, 2, 6, 5]], axis=0)
            ax.plot(
                [center_bottom[0], center_bottom_forward[0]],
                [center_bottom[1], center_bottom_forward[1]],
                c="y",
                lw=3,
            )

        for pts_score in keypoints:
            pts = pts_score[:8]
            if divide:
                pts = pts / 8.0
            for i in range(4):
                ax.plot(pts[2 * i + 1], pts[2 * i + 0], "r*")
            ax.plot([pts[1], pts[3]], [pts[0], pts[2]], c="y", lw=5)
            ax.plot([pts[3], pts[5]], [pts[2], pts[4]], c="g", lw=5)
            ax.plot([pts[5], pts[7]], [pts[4], pts[6]], c="b", lw=5)
            ax.plot([pts[7], pts[1]], [pts[6], pts[0]], c="r", lw=5)

            # ax.plot([pts[7], pts[7]+10*pts_score[11]], [pts[6], pts[6]+10*pts_score[12]], c='c', lw=2)
            # annotations = pts_score[13]
            # ax.annotate(str(annotations), xy=(pts[7]+10*pts_score[11], pts[6]+10*pts_score[12]), xytext=(pts[7]+10*pts_score[11], pts[6]+10*pts_score[12]), arrowprops=dict(facecolor='black', shrink=0.05),)

        ax.axis("off")
    if save_filename:
        plt.savefig(save_filename)
    else:
        plt.show()
    plt.close()