def plotTwoContours(contour, contour2, diffmap, mapName1="map1", mapName2="map2", colormap="jet", plot_scalebar=False):
        #DEBUG:
        #max_range = 0.07
        max_range = max(diffmap.max(),abs(diffmap.min()))
	norm = colors.Normalize(vmin=-max_range, vmax=max_range) 
	print('max_range is '+str(max_range))
	title = mapName1+' - '+mapName2
	filename = title+'.pdf'
	fig = plt.figure()
	plt.title(title)
        plt.hold(True)
	plt.imshow(diffmap, origin='upper', norm=norm)
        ax = plt.axes()
        if plot_scalebar:
            addScaleBar(ax)
        else:
            scaleTicks(ax)
        #if(colormap == "rwb")
        plt.set_cmap(colormap)
	plt.colorbar()
	plt.contour(contour, [0], colors='r')
	#plt.contour(image2, [0], colors='b')
	plt.contour(contour2, [0], colors='b')
	plt.hold(False)
	plt.savefig(filename)
        #plt.show()
	return diffmap

def main():
	for unitFile in os.listdir(sourceFolder):
		if os.path.isdir(sourceFolder+unitFile):			
			unitName = unitFile.rsplit('_', 1)[0]
			print unitName
			staMatrixFile = scipy.io.loadmat(sourceFolder+unitFile+'/stavisual_lin_array_'+unitName+'.mat')
			staMatrix = staMatrixFile['STAarray_lin']
			xLength = staMatrix.shape[0]
			yLength = staMatrix.shape[1]
			zLength = staMatrix.shape[2]
			for zAxis in range(zLength):
				print 'desde disco'
				
				fig = plt.figure()
				fig.set_size_inches(1, 1)
				data = staMatrix[:,:,zAxis]
				#plt.pcolormesh( staMatrix[:,:,zAxis],vmin = 0,vmax = 255, cmap=cm.jet )
				ax = plt.Axes(fig, [0., 0., 1., 1.])
				ax.set_axis_off()
				fig.add_axes(ax)
				plt.set_cmap(cm.jet)
				ax.imshow(data, aspect = 'auto')
				plt.savefig(outputFolder+unitName+str(zAxis)+".png",format='png',dpi=31)
				plt.close()
			
	return 0

def dim_sensitivity_plot(x, Y, fname, show_legend=True):

    plt.rc('text', usetex=True)
    plt.rc('font', family='serif')

    plt.figure(figsize=(3, 3))
    plt.xlabel('$d$', size=FONTSIZE)
    plt.ylabel('ROC AUC', size=FONTSIZE)

    plt.set_cmap('Set2')

    lines = []
    for i, label in enumerate(KEYS):
        line_data = Y.get(label)

        if line_data is None:
            continue
        
        line, = plt.plot(x, line_data, label=label, marker=MARKERS[i],
                         markersize=0.5 * FONTSIZE, color=COLORS[i])
        lines.append(line)



    if show_legend:
        plt.legend(handles=lines)
        plt.legend(loc='lower right')
    plt.xscale('log', basex=2)
    plt.xticks(x, [str(y) for y in x], size=FONTSIZE)
    plt.yticks(size=FONTSIZE)
    plt.tight_layout()

    plt.savefig(fname)

    def plot_summary(self):
        import matplotlib
        import matplotlib.pyplot as plt

        fig = plt.figure(figsize=(10,8)) 
        ax = plt.axes((0.08, 0.08, 0.87, 0.80)) 
        plt.set_cmap("Spectral") 
        
        # Get all candidates and sort by sigma
        allcands = self.get_all_cands()
        sigmas = Num.array([c.sigma for c in allcands])
        isort = sigmas.argsort()
        sigmas = sigmas[isort]
        freqs = Num.array([c.f for c in allcands])[isort]
        dms = Num.array([c.DM for c in allcands])[isort]
        numharms = Num.array([c.numharm for c in allcands])[isort]

        # Plot the all candidates 
        plt.scatter(freqs, dms, s=8+sigmas**1.7, c=Num.log2(numharms), \
                                marker='o', alpha=0.7, zorder=-1) 
  
        # Add colorbar 
        fmtr = matplotlib.ticker.FuncFormatter(lambda x, pos: "%d" % 2**x) 
        cb = plt.colorbar(ticks=(0,1,2,3,4), format=fmtr) 
        cb.set_label("Num harmonics summed") 
         
        plt.xscale('log', base=10.0) 
        plt.xlim(0.5, 10000) 
        plt.ylim(-10, 1200) 
        plt.xlabel("Freq (Hz)") 
        plt.ylabel(r"DM (pc cm$^{-3}$)") 
        return fig

def _vis_readout_maps(outputs, global_step, output_dir, metric_summary, N):
  # outputs is [gt_map, pred_map]:
  if N >= 0:
    outputs = outputs[:N]
  N = len(outputs)

  plt.set_cmap('jet')
  fig, axes = utils.subplot(plt, (N, outputs[0][0].shape[4]*2), (5,5))
  axes = axes.ravel()[::-1].tolist()
  for i in range(N):
    gt_map, pred_map = outputs[i]
    for j in [0]:
      for k in range(gt_map.shape[4]):
        # Display something like the midpoint of the trajectory.
        id = np.int(gt_map.shape[1]/2)

        ax = axes.pop();
        ax.imshow(gt_map[j,id,:,:,k], origin='lower', interpolation='none',
                  vmin=0., vmax=1.)
        ax.set_axis_off();
        if i == 0: ax.set_title('gt_map')

        ax = axes.pop();
        ax.imshow(pred_map[j,id,:,:,k], origin='lower', interpolation='none',
                  vmin=0., vmax=1.)
        ax.set_axis_off();
        if i == 0: ax.set_title('pred_map')

  file_name = os.path.join(output_dir, 'readout_map_{:d}.png'.format(global_step))
  with fu.fopen(file_name, 'w') as f:
    fig.savefig(f, bbox_inches='tight', transparent=True, pad_inches=0)
  plt.close(fig)

def getHash(im_file):
    im = np.array(Image.open(im_file).convert('L'))
    f=np.fft.fft2(im)
    rmax,cmax=f.shape    
    sg_r=np.zeros((2*wd,wd))
    sg_r[0:wd,:]=np.abs(f.real[0:wd,0:wd])
    sg_r[wd:2*wd,:]=np.abs(f.real[rmax-wd:rmax,0:wd])
    sg_i=np.zeros((2*wd,wd))
    sg_i[0:wd,:]=np.abs(f.imag[0:wd,0:wd])
    sg_i[wd:2*wd,:]=np.abs(f.imag[rmax-wd:rmax,0:wd])
    
    bsg=np.zeros((2*wd,2*wd-2),dtype=bool)
    bsg[:,0:wd-1]=sg_r[:,0:wd-1]<sg_r[:,1:wd]
    bsg[:,wd-1:2*wd-2]=sg_i[:,0:wd-1]<sg_i[:,1:wd]
    
    return bsg
    
    f2=np.zeros((rmax,cmax))+1j*np.zeros((rmax,cmax))
    
    f2.real[0:wd,0:wd-1]=bsg[0:wd,0:wd-1]
    f2.imag[wd:2*wd,0:wd-1]=bsg[0:wd,wd-1:2*wd-2]
    f2.real[rmax-wd:rmax,0:wd-1]=bsg[wd:2*wd,0:wd-1]
    f2.imag[rmax-wd:rmax,0:wd-1]=bsg[wd:2*wd,wd-1:2*wd-2]
    #plt.close("all")
    plt.figure()
    plt.imshow(bsg,interpolation='none')
    plt.set_cmap("gray")

def produce_regions(masks, visualize=False):
    """given the proposal segmentation masks for an image as a [width, height, proposal_num]
    matrix outputs all regions in the image"""
    width, height, n_prop = masks.shape
    t = ('u8,'*int(np.math.ceil(float(n_prop) / 64)))[:-1]
    bv = np.zeros((width, height), dtype=np.dtype(t))
        
    for i in range(n_prop):
        m = masks[:, :, i]
        a = 'f%d' % (i / 64)    
        h = m * np.long(2 ** (i % 64))
        if n_prop >= 64:
            bv[a] += h
        else:
            bv += h


    un = np.unique(bv)
    regions = np.zeros((width, height), dtype="uint16")
    for i, e in enumerate(un):
        regions[bv == e] = i
    if visualize:
        plt.figure()
        plt.imshow(regions)
        plt.set_cmap('prism')
        plt.colorbar()
    return regions

    def plot_reals(self, nr=25, hardcopy=0, hardcopy_filename='reals', nanval=-997799, filternan=1):
        import matplotlib.pyplot as plt
        import matplotlib.gridspec as gridspec
        import numpy as np
        from scipy import squeeze
        
        
        nx, ny, nz = self.par['simulation_grid_size']
        
        nr = np.min((self.par['n_real'], nr))
        nsp = int(np.ceil(np.sqrt(nr)))

        fig = plt.figure(1)
        sp = gridspec.GridSpec(nsp, nsp, wspace=0.1, hspace=0.1)
        plt.set_cmap('hot')
        for i in range(0, nr):
            ax1 = plt.Subplot(fig, sp[i])
            fig.add_subplot(ax1)
            if filternan==1:
                self.sim[i][self.sim[i]==nanval] = np.nan
                
            D=squeeze(np.transpose(self.sim[i]))
            plt.imshow(D, extent=[self.x[0], self.x[-1], self.y[0], self.y[-1]], interpolation='none')
            plt.title("Real %d" % (i + 1))

        fig.suptitle(self.method + ' - ' + self.parameter_filename, fontsize=16)
        if (hardcopy):
            plt.savefig(hardcopy_filename)

        plt.show(block=False)

    def plot1(self, img, path=None, mask=None, interpolation='none', colorMap='jet', suffix=''):
        import matplotlib.pyplot as plt
        if path != None:
            plt.ioff()

        if isinstance(mask, np.ndarray):
            img = img[:,:] * mask

        plt.imshow(img, interpolation=interpolation)
        plt.set_cmap(colorMap)
        cbar = plt.colorbar()
        cbar.set_ticks([])
        if path != None:
            if suffix == None:
                fout = osp.join(path, '{0}.png'.format(self.label))
            else:
                fout = osp.join(path, '{0}_{1}.png'.format(self.label, suffix))
            try:
                plt.savefig(fout)
            except IOError:
                raise IOError('in classifiers.output, no such file or directory: {0}'.format(path))
        else:
            if suffix == None:
                plt.title('{0}'.format(self.label))
            else:
                plt.title('{0} - {1}'.format(self.label, suffix))
            plt.show()

        plt.close()

 def set_cmap(self, i):
     """Set colormap."""
     if 0 <= i < len(self.cmaps):
         name = self.cmaps[i]
         if self.is_reverse_cmap:
             name = dwi.plot.reverse_cmap(name)
         plt.set_cmap(name)

def tracestackedslab(infile,depthinc=5,llinc=((6371*np.pi)/360),cval=0.5):
	'''
	Takes a netcdf file containing a stacked, normed profiles and attempts to contour what might be a slab - can use this to estimate dip, profile etc
	The values of depthinc and llinc are defaults from the Ritsema code, which creates slices over angles of 180 degrees
	and with a depth increment of 5 km

	This produces a map showing the slice and contours in stacked velocity perturbation at a chosen level (typically 0.5)
	'''

	Mantlebase = 2895

	infile = Dataset(infile, model='r')

	filevariables = infile.variables.keys()

	#Get data from the netCDF file
	depths = infile.variables['y'][:]
	lengths = infile.variables['x'][:]
	data = infile.variables['z'][:][:]

	infile.close()

	#print np.shape(data)
	#print np.shape(lengths)
	#print np.shape(depths)

	#Use image processing suite to find contours
	contours = measure.find_contours(data,cval)

	#Various plotting commands to produce the figure
	fig, ax = plt.subplots()

	thousandkm = int((Mantlebase-1000)/depthinc)
	sixsixtykm = int((Mantlebase-660)/depthinc)

	plt.set_cmap('jet_r')

	ax.imshow(data, interpolation='linear',aspect='auto')

	ax.plot([0,len(lengths)],[thousandkm,thousandkm],'k--',label='1000km')
	ax.plot([0,len(lengths)],[sixsixtykm,sixsixtykm],'k-',label='660km')

	for n, contour in enumerate(contours):
		if n == 0:
			ax.plot(contour[:, 1], contour[:, 0], 'r-', linewidth=2,label='contour at %g' %cval)
		else:
			ax.plot(contour[:, 1], contour[:, 0], 'r-', linewidth=2)


	ax.set_ylim([0,len(depths)])
	ax.set_xlim([len(lengths),0])
	ax.set_title('Stacked slab image from netCDF')
	plt.xlabel('Cross section x increment')
	plt.ylabel('Cross section depth increment')
	plt.legend(loc='best')

	#plt.gca().invert_yaxis()
	#plt.gca().invert_xaxis()

	plt.show(block=False)

def testModel(modelFilename, testData, size):
    

    # Create shared data
    test_x, test_y = make_shared(testData, size)
    
    # Load the model
    persistence = PersistenceManager()
    persistence.set_filename(modelFilename)
    x, y, classifier = persistence.load_model()
    
    # Create prediction function.
    pred_func = theano.function(inputs = [],
                                outputs = [classifier.y_pred, classifier.errors(y)],
                                givens = {x:test_x, y:test_y})    


    preds = pred_func()
    print("prediction done")

    plt.set_cmap('gray')    
    for i in range(len(testData)):
        s = testData[i]
        print("Dist: {0}".format(preds[1][i]))        
        normalizedCoords = preds[0][i]
        px = np.round(normalizedCoords[0] * s.width)
        py = np.round(normalizedCoords[1] * s.height)

        pxCopy = s.getAnnotated()
        pxCopy[py, px] = 0.85*np.max(pxCopy)

        plt.imshow(pxCopy)
        plt.show()    

def visualize_matrix(matrix, n_imgs, imgsize, outputfile, cmap="gray", dpi=150):
    """

    """
    n_h,n_w = layout_shape(n_imgs)
    
    receptive_fields = numpy.zeros((n_h * imgsize, n_w * imgsize), dtype=theano.config.floatX)
    
    for i in range(n_h):
        for j in range(n_w):
            img = matrix[i * n_w + j]
            img = img.reshape((imgsize, imgsize))
            receptive_fields[i * imgsize: (i + 1) * imgsize, j * imgsize: (j + 1) * imgsize] = img
    
    fig = plt.figure()
    
    ax = fig.add_subplot(1,1,1)
    xticks = numpy.arange(0, n_w * imgsize, imgsize)                                              
    yticks = numpy.arange(0, n_h * imgsize, imgsize)
    ax.set_xticks(xticks)
    ax.set_xticklabels([i for (i,x) in enumerate(xticks)])
    ax.set_yticks(yticks)
    ax.set_yticklabels([i for (i,y) in enumerate(yticks)])
    ax.grid(which="both", linestyle='-')
    
    plt.set_cmap(cmap)
    plt.imshow(receptive_fields)
    plt.savefig(outputfile, dpi=dpi)

def plot_topography(grid, elev):
    ''' 
    This function takes the DEM read in below and plots the
    elevations for visualization purposes.
    '''
    # Get a 2D array version of the elevations for plotting purposes
    elev_raster = grid.node_vector_to_raster(elev,True)
    
    # Everything below plots the topography and sampling points
    levels = []
    # To better create a colorbar...
    x_up = 2475
    while x_up !=2600:
        levels.append(x_up)
        x_up+=1
    plt.figure('Topography')
    plt.contourf(elev_raster, levels, colors='k')
    plt.set_cmap('bone')
    plt.colorbar()
    
    # To plot the study node and outlet node on the DEM...
    plt.plot([122],[140],'cs', label= 'Study Node')
    plt.plot([152],[230], 'wo', label= 'Outlet')
    
    plt.legend(loc=3)

def test_mask_loss_sobel():
    th_mask, th_img = T.tensor4(), T.tensor4()
    ml = mask_loss_sobel(th_mask, th_img)
    mask_loss = theano.function([th_mask, th_img],
                                [ml.loss] + list(ml.sobel_mask) +
                                list(ml.sobel_img))

    mask_idx = next(masks(1))
    image_ok = 0.5 * np.ones_like(mask_idx)
    image_ok[mask_idx > MASK["IGNORE"]] = 1
    image_ok[mask_idx < MASK["BACKGROUND_RING"]] = 0

    print()
    loss, sobel_mask_x, sobel_mask_y, sobel_img_x, sobel_img_y = \
        mask_loss(mask_idx, image_ok)
    plt.set_cmap('gray')
    plt.subplot(221)
    plt.imshow(sobel_mask_x[0, 0])
    plt.subplot(222)
    plt.imshow(sobel_mask_y[0, 0])
    plt.colorbar()
    plt.subplot(223)
    plt.imshow(sobel_img_x[0, 0])
    plt.subplot(224)
    plt.imshow(sobel_img_y[0, 0])
    plt.colorbar()
    plt.savefig("mask_loss_sobel.png")
    print()
    print("mask_loss: {}".format(mask_loss(mask_idx, image_ok)))
    assert loss == 0

 def _plot_single_map1(self, path, cmap, signo, dist_map, threshold, constrained, stretch, colorMap, suffix):
     import matplotlib.pyplot as plt
     if path != None:
         plt.ioff()
     grad = self.get_single_map(signo, cmap, dist_map, threshold, constrained, stretch)
     plt.imshow(grad, interpolation='none')
     plt.set_cmap(colorMap)
     cbar = plt.colorbar()
     cbar.set_ticks([])
     if path != None:
         if suffix == None:
             fout = osp.join(path, '{0}_{1}.png'.format(self.label, signo))
         else:
             fout = osp.join(path, '{0}_{1}_{2}.png'.format(self.label, signo, suffix))
         try:
             plt.savefig(fout)
         except IOError:
             raise IOError('in classifiers.output, no such file or directory: {0}'.format(path))
     else:
         if suffix == None:
             plt.title('{0} - EM{1}'.format(self.label, signo))
         else:
             plt.title('{0} - EM{1} - {2}'.format(self.label, signo, suffix))
         plt.show()
     plt.close()

def show_network(images, side_length, sample=1):
    """
    Expects images to be (n_pixels, n_images).  n_pixels should be == side_length ** 2
    """
    n_images = int(images.shape[1] * sample)
    cols = int(math.sqrt(n_images))
    rows = math.ceil(n_images / cols)
    # padding
    image_size = side_length + 1
    output = np.zeros((image_size * rows, image_size * cols))
    image_mask = np.random.randint(0, images.shape[1], n_images)
    norm = Normalize()
    for i, img in enumerate(image_mask):
        this_image = images[:, img].reshape((side_length, side_length)).copy()
        # Center and normalize
        this_image -= this_image.mean()
        this_image = norm(this_image)
        # Get offsets
        offset_col = image_size * (i % cols)
        offset_col_end = offset_col + image_size - 1
        offset_row = image_size * (math.floor(i / cols))
        offset_row_end = offset_row + image_size - 1
        output[offset_row:offset_row_end, offset_col:offset_col_end] = this_image

    pyplot.imshow(output)
    pyplot.set_cmap(cm.gray)
    pyplot.show()

 def show(self):
     import matplotlib.pyplot as plt
     print("number of images: {}".format(len(self.I)))    
     for i in xrange(len(self.jsonfiles)):
         f=open(self.jsonfiles[i],"r")
         js=json.loads(f.read())
         f.close()
         
         ##init and show
         img_path=''
         if js['path'][0:2]=='./':
             img_path= rootdir + js['path'][1:]
         elif js['path'][0]=='/':
             img_path= rootdir + js['path']
         else:
             img_path= rootdir + '/' +js['path'] 
             
         print(img_path)
         im=np.array(Image.open(img_path).convert('L'))
         plt.hold(False)        
         plt.imshow(im)
         plt.hold(True)
         for j in range(self.size):
             #samples[x]=[0_class,1_img, 2_row, 3_column]^T
             if self.samples[1,j]==i:
                 plt.text(self.samples[3,j], self.samples[2,j], "%03d"%self.samples[0,j], fontsize=12,color='red')
                 #plt.plot(self.samples[3,j],self.samples[2,j],markers[self.samples[0,j]])
         plt.set_cmap('gray')
         plt.show()
         plt.ginput(1)
     plt.close('all')

def compare_spectra():
    import mywfc3.stgrism as st
    import unicorn
    
    ### Fancy colors
    import seaborn as sns
    import matplotlib.pyplot as plt
    cmap = sns.cubehelix_palette(as_cmap=True, light=0.95, start=0.5, hue=0.4, rot=-0.7, reverse=True)
    cmap.name = 'sns_rot'
    plt.register_cmap(cmap=cmap)
    sns.set_style("ticks", {"ytick.major.size":3, "xtick.major.size":3})
    plt.set_cmap('sns_rot')
    #plt.gray()
    
    fig = st.compare_methods(x0=787, y0=712, v=np.array([-1.5,4])*0.6, NX=180, NY=40, direct_off=100, final=True, mask_lim = 0.02)
    #fig.tight_layout()
    unicorn.plotting.savefig(fig, '/tmp/compare_model_star.pdf', dpi=300)

    fig = st.compare_methods(x0=485, y0=332, v=np.array([-1.5,4])*0.2, NX=180, NY=40, direct_off=100, final=True, mask_lim = 0.1)
    unicorn.plotting.savefig(fig, '/tmp/compare_model_galaxy.pdf', dpi=300)

    fig = st.compare_methods(x0=286, y0=408, v=np.array([-1.5,4])*0.08, NX=180, NY=40, direct_off=100, final=True, mask_lim = 0.1)
    unicorn.plotting.savefig(fig, '/tmp/compare_model_galaxy2.pdf', dpi=300)

    fig = st.compare_methods(x0=922, y0=564, v=np.array([-1.5,4])*0.2, NX=180, NY=40, direct_off=100, final=True, mask_lim = 0.15)
    unicorn.plotting.savefig(fig, '/tmp/compare_model_galaxy3.pdf', dpi=300)

def displayData(X):
    width = 20
    rows, cols = 10, 10
    out = zeros((width * rows, width * cols))
    m = X.shape[0]  # データ数

    # 5000個のデータセットから適当に100個選ぶ
    rand_indices = random.permutation(m)[0: rows * cols]

    counter = 0
    for y in range(0, rows):
        for x in range(0, cols):
            start_x = x * width
            start_y = y * width
            out[start_x: start_x + width, start_y: start_y + width] = X[rand_indices[counter]].reshape(width, width).T
            counter += 1

    img = scipy.misc.toimage(out)
    figure = pyplot.figure()
    pyplot.tick_params(labelbottom="off")
    pyplot.tick_params(labelleft="off")
    pyplot.set_cmap(pyplot.gray())
    axes = figure.add_subplot(111)
    axes.imshow(img)
    pyplot.savefig("digits.png")