def color_plot(targets,k_DD,f_DD,coherence,azimuths):

    f_DD_boundaries = ( f_DD - (f_DD[1]-f_DD[0])/2 )
    d_phi = 1/coherence.shape[1]
    angles_over_pi = np.arange(0,1+d_phi,1/coherence.shape[1])

    plt.pcolor(angles_over_pi,f_DD_boundaries,coherence)
    if show_azimuths:
        for azimuth in azimuths:
            plt.axvline(mod(azimuth),color="k",linewidth=2)

    plt.xlim(0,1)
    plt.ylim(0,f_DD_boundaries[-1])
    plt.clim(-1,1)

    plt.xlabel(r"$\phi_{DD}/\pi$")
    plt.ylabel("$f_{}$".format(k_DD))

    cbar = plt.colorbar()
    cbar.set_label("Coherence")

    plt.tight_layout()
    plt.savefig(figname(targets))

    if show: plt.show()

	def on_keypress(self,event):
		global colmax
		global colmin
		if event.key in ['1', '2', '3', '4', '5', '6', '7','8', '9', '0']:
			if not os.path.exists(write_dir + runtag):
				os.mkdir(write_dir + runtag)
			recordtag = write_dir + runtag + "/" + runtag + "_" + event.key + ".txt"
			print "recording filename in " + recordtag
			f = open(recordtag, 'a+')
			f.write(self.filename+"\n")
			f.close()
		if event.key == 'p':
			if not os.path.exists(write_dir + runtag):
				os.mkdir(write_dir + runtag)
			pngtag = write_dir + runtag + "/%s.png" % (self.filename)	
			print "saving image as " + pngtag 
			P.savefig(pngtag)
		if event.key == 'e':
			if not os.path.exists(write_dir + runtag):
				os.mkdir(write_dir + runtag)
			epstag = write_dir + runtag + "/%s.eps" % (self.filename)	
			print "saving image as " + epstag 
			P.savefig(epstag, format='eps')
		if event.key == 'r':
			colmin = self.inarr.min()
			colmax = self.inarr.max()
			P.clim(colmin, colmax)
			P.draw()

def showMechStiff(model, coords, *args, **kwargs):
    """Show mechanical stiffness matrix using :func:`~matplotlib.pyplot.imshow`.
    By default, *origin=lower* keyword  arguments are passed to this function, 
    but user can overwrite these parameters."""

    import math
    import matplotlib
    import matplotlib.pyplot as plt
    arange = np.arange(model.numAtoms())
    model.buildMechStiff(coords)

    if not 'origin' in kwargs:
        kwargs['origin'] = 'lower'
    if 'jet_r' in kwargs:
        import matplotlib.cm as plt
        kwargs['jet_r'] = 'cmap=cm.jet_r'
        
    MechStiff = model.getStiffness()
    matplotlib.rcParams['font.size'] = '14'
    fig = plt.figure(num=None, figsize=(10,8), dpi=100, facecolor='w')
    show = plt.imshow(MechStiff, *args, **kwargs), plt.colorbar()
    plt.clim(math.floor(np.min(MechStiff[np.nonzero(MechStiff)])), \
                                           round(np.amax(MechStiff),1))
    #plt.title('Mechanical Stiffness Matrix')# for {0}'.format(str(model)))
    plt.xlabel('Indices', fontsize='16')
    plt.ylabel('Indices', fontsize='16')
    if SETTINGS['auto_show']:
        showFigure()
    return show

def my_plot(c_plot, step):
    filename = "chs-step" + str(step) + ".png"
    plt.imshow(c_plot)
    plt.colorbar()
    plt.clim(0, 1)
    plt.savefig(filename)
    plt.clf()

def create_animation_rate(rate, frames=100, interval=1, fps=10, dpi=90, filename='default.mp4', format='.mp4'):
    """
    Documentation needed 
    """
    
    filename = filename + '-rate_animation' + format
    # Initiate figure 
    fig = plt.figure()
    ims = plt.imshow(rate[::,::,1])
    plt.xlabel('Neuron\'s x coordinate')
    plt.ylabel('Neuron\'s y coordinate')
    plt.title('Firing rate in a 20 ms window')

    cbar = plt.colorbar()
    cbar.ax.set_ylabel('Firing Rate (Hz)')
    plt.clim(0,70)
    
    # Define how to update frames 
    def updatefig(i):
        ims.set_array( rate[:,:,i] )
        return ims,

    # run and save the animation
    image_animation = animation.FuncAnimation(fig,updatefig, frames=frames, interval=interval, blit = True)
    image_animation.save(filename, fps=fps, dpi=dpi)
    plt.show()

def _colormap_plot_array_response(cmaps):
    """
    Plot for illustrating colormaps: array response.

    :param cmaps: list of :class:`~matplotlib.colors.Colormap`
    :rtype: None
    """
    import matplotlib.pyplot as plt
    from obspy.signal.array_analysis import array_transff_wavenumber
    # generate array coordinates
    coords = np.array([[10., 60., 0.], [200., 50., 0.], [-120., 170., 0.],
                       [-100., -150., 0.], [30., -220., 0.]])
    # coordinates in km
    coords /= 1000.
    # set limits for wavenumber differences to analyze
    klim = 40.
    kxmin = -klim
    kxmax = klim
    kymin = -klim
    kymax = klim
    kstep = klim / 100.
    # compute transfer function as a function of wavenumber difference
    transff = array_transff_wavenumber(coords, klim, kstep, coordsys='xy')
    # plot
    for cmap in cmaps:
        plt.figure()
        plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
                   np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
                   transff.T, cmap=cmap)
        plt.colorbar()
        plt.clim(vmin=0., vmax=1.)
        plt.xlim(kxmin, kxmax)
        plt.ylim(kymin, kymax)
    plt.show()

def main(unused_argv):
    global heat_map
    tf_records_file_names = sorted(os.listdir(utils.HEAT_MAP_TF_RECORDS_DIR))
    print(tf_records_file_names)
    tf_records_file_names = tf_records_file_names[2:3]
    for wsi_filename in tf_records_file_names:
        print('Generating heatmap for: %s' % wsi_filename)
        tf_records_dir = os.path.join(utils.HEAT_MAP_TF_RECORDS_DIR, wsi_filename)
        raw_patches_dir = os.path.join(utils.HEAT_MAP_RAW_PATCHES_DIR, wsi_filename)
        heatmap_rgb_path = os.path.join(utils.HEAT_MAP_WSIs_PATH, wsi_filename)
        assert os.path.exists(heatmap_rgb_path), 'heatmap rgb image %s does not exist' % heatmap_rgb_path
        heatmap_rgb = Image.open(heatmap_rgb_path)
        heatmap_rgb = np.array(heatmap_rgb)
        heatmap_rgb = heatmap_rgb[:, :, :1]
        heatmap_rgb = np.reshape(heatmap_rgb, (heatmap_rgb.shape[0], heatmap_rgb.shape[1]))
        heat_map = np.zeros((heatmap_rgb.shape[0], heatmap_rgb.shape[1]), dtype=np.float32)
        assert os.path.exists(raw_patches_dir), 'raw patches directory %s does not exist' % raw_patches_dir
        num_patches = len(os.listdir(raw_patches_dir))
        assert os.path.exists(tf_records_dir), 'tf-records directory %s does not exist' % tf_records_dir
        dataset = Dataset(DATA_SET_NAME, utils.data_subset[4], tf_records_dir=tf_records_dir, num_patches=num_patches)
        build_heatmap(dataset)
        # Image.fromarray(heat_map).save(os.path.join(utils.HEAT_MAP_DIR, wsi_filename), 'PNG')
        plt.imshow(heat_map, cmap='hot', interpolation='nearest')
        plt.colorbar()
        plt.clim(0.00, 1.00)
        plt.axis([0, heatmap_rgb.shape[1], 0, heatmap_rgb.shape[0]])
        plt.savefig(str(os.path.join(utils.HEAT_MAP_DIR, wsi_filename))+'_heatmap.png')
        plt.show()

def plotSVM(solver, lowerBound, upperBound, step):
    
    assert lowerBound < upperBound

    X = arange(lowerBound, upperBound, step)
    Y = arange(lowerBound, upperBound, step)
    X,Y = meshgrid(X,Y)
    Z = zeros(X.shape)

    for i in range(len(X)):
        for j in range(len(Y)):
            
            #Classify the result
            result = int( svm_predict([0.], [[ X[i][j], Y[i][j] ]], solver, '-q')[0][0] )
            if result == 0:
                Z[i][j] = 0 #lower limit
            else:
                Z[i][j] = 100 #higher limit

    plot.imshow(Z)
    plot.gcf()
    plot.clim()
    plot.title("SVM Activation")

    return plot

def main():
    import matplotlib.pyplot as plt
    x, y = np.mgrid[-1:1:200j, -1:1:200j]
    rho = np.sqrt(x**2 + y**2)
    theta = np.arctan(y / x)
    theta[(theta.shape[0] / 2):] += np.pi
    nm_pairs = list(n_m_upto(4, 2))
    plt.ioff()
    cur_max = -np.inf
    cur_min = np.inf
    main_axes = plt.gca()
    for index, (n, m) in enumerate(nm_pairs):
        image = zernike_exp(rho, theta, n, m)
        cur_max = max(cur_max, np.max(np.real(image[~np.isnan(image)])),
                      np.max(np.imag(image[~np.isnan(image)])))
        cur_min = min(cur_min, np.min(np.real(image[~np.isnan(image)])),
                      np.min(np.imag(image[~np.isnan(image)])))
        plt.subplot(2, len(nm_pairs), index + 1)
        plt.imshow(np.real(image), cmap=plt.cm.gray)
        plt.title('$\\mathrm{real}(V_{%d,%d}(\\rho, \\theta))$' % (n, m))
        plt.axis('off')
        plt.subplot(2, len(nm_pairs), len(nm_pairs) + index + 1)
        plt.imshow(np.imag(image), cmap=plt.cm.gray)
        plt.title('$\\mathrm{imag}(V_{%d,%d}(\\rho, \\theta))$' % (n, m))
        plt.axis('off')
    for index in range(2 * len(nm_pairs)):
        plt.subplot(2, len(nm_pairs), index + 1)
        plt.clim(cur_min, cur_max)
    print "cur_min =", cur_min
    print "cur_max =", cur_max
    plt.show()
    plt.ion()

def plot_fld(fldname, step, cra=None, xlim=None, ylim=None):
    flds = psc.PscFields(path, step, "p")
    fld = flds[fldname]
    crd_nc = flds.crd_nc

    plt.figure()
    kwargs = {}

    plt.pcolormesh(crd_nc[1], crd_nc[2], fld[:,:,0], **kwargs)
    if cra:
        plt.clim(*cra)
    plt.colorbar()
    plt.title(fldname)
    plt.xlabel("y")
    plt.ylabel("z")

    if xlim:
        plt.xlim(*xlim)
    else:
        plt.xlim(crd_nc[1][0], crd_nc[1][-1])

    if ylim:
        plt.ylim(*ylim)
    else:
        plt.ylim(crd_nc[2][0], crd_nc[2][-1])

    plt.axes().set_aspect('equal')
    plt.savefig("%s-%06d.png" % (fldname, step), dpi=100)
    plt.close()

def data_movie(data, fig=None, pause=None, **kwargs):
    """
    Display all images of a data cube as a movie.

    Arguments
    ---------
    data: ndarray
      A data cube. Third dimension should be the image index.
    fig: Figure instance (optional).
       A figure instance where the images will be displayed.
    kwargs: Keywor arguments.
       Other keyword arguments are passed to imshow.

    Returns
    -------
    Nothing.
    """
    if fig is None:
        fig = plt.figure()
    a = fig.gca()
    dmin, dmax = data.min(), data.max()
    n = data.shape[-1]
    im0 = a.imshow(data[..., 0].T, origin="lower", **kwargs)
    plt.draw()
    for k in xrange(n):
        im0.set_data(data[..., k].T)
        plt.title("Image " + str(k))
        plt.clim([dmin, dmax])
        plt.draw()
        if pause is not None:
            time.sleep(pause)

 def plot_grammar_matrices(self,
                           folder_path,
                           folder_name,
                           A_matrix = np.zeros(0),
                           B_matrix = np.zeros(0)):
     if folder_name in os.listdir(folder_path):
         shutil.rmtree(folder_path + '/' + folder_name,
                       True)
     os.mkdir(folder_path + '/' + folder_name)
     if(len(A_matrix) == 0):
         A_matrix = self.A
     if(len(B_matrix) == 0):
         B_matrix = self.B
     assert(A_matrix.shape[0] == A_matrix.shape[1] == A_matrix.shape[2] == B_matrix.shape[0])
     N = A_matrix.shape[0]
     for i in xrange(N):
         plt.subplot(211)
         plt.title('A %d' % i)
         plt.imshow(A_matrix[i], interpolation = 'None')
         plt.clim(0, 1.0)
         plt.subplot(212)
         plt.plot(range(len(B_matrix[i])), B_matrix[i], linestyle = 'None', marker = 'o')
         plt.ylim(-0.2, 1.0)
         plt.xlim(-1, len(B_matrix[i]))
         plt.title('B %d' % i)
         plt.savefig(folder_path + '/' + folder_name + '/' + string.lower(folder_name) + '_rule_' + str(i) + '.png', dpi = 300)
         plt.close()

def animate(i):
    plt.clf()
    m, _ = atm.pcolor_latlon(animdata[i], axlims=axlims, cmap=cmap)
    plt.clim(cmin, cmax)
    day = animdata[daynm + 'rel'].values[i]
    plt.title('%s %s RelDay %d' % (var, yearstr, day))
    return m

 def heat_map(self,DL):
   fig = plt.figure()
   ax = fig.add_subplot(1,1,1)
   plt.imshow(DL)
   plt.clim(DL.mean()-3*DL.std(),DL.mean()+3*DL.std())
   plt.colorbar()
   plt.show(block=False)

def plot_gridsearch_result(tau_sses, tau_ranges, tau_opt, ptitle='default title',save=False):
    fig = plt.figure(figsize=(24,20), facecolor='white')
    subrows = np.ceil(np.sqrt(tau_ranges[-1].shape))
    subcols = np.ceil(tau_ranges[-1].shape/subrows)
    gs = gridspec.GridSpec(int(subrows), int(subcols), left=0.1, right=0.90, bottom=0.1, top=0.9)
    gs.update(hspace=0.5, wspace=0.5)
    
    for ctau3_idx,ctau3 in enumerate(tau_ranges[-1]):
        plt.subplot(gs[ctau3_idx])
        img = plt.imshow(tau_sses[:,:,ctau3_idx], aspect='auto', interpolation='none', origin='lower')
        img.axes.set_yticks(np.arange(len(tau_ranges[0])))
        img.axes.set_yticklabels(tau_ranges[0])
        img.axes.set_xticks(np.arange(len(tau_ranges[1]),step=2))
        img.axes.set_xticklabels(np.arange(tau_ranges[1][0],np.max(tau_ranges[1]),step=10))
        plt.xlabel('Second time constant [ms]',fontsize=12)
        plt.ylabel('First time constant [ms]',fontsize=12)
        plt.title('Total SSE, for tau 3 = '+str(ctau3) + 'ms',fontsize=12)
        plt.colorbar()
        plt.clim(tau_sses.min(),tau_sses.min()*1.5)
        
    plt.suptitle(ptitle,fontsize=20)
    
    fig.text(0.8,0.95, 'Optimal set of taus: '+str(tau_opt.squeeze()))
    
    if save==True:
        if not os.path.exists(figure_path):
            os.makedirs(figure_path)
            
        plt.savefig(figure_path + ptitle, dpi=120)
        plt.close(fig)

 def plot(self,vmin=np.NaN,vmax=np.NaN,colorbarcenter=False,colorbar=cm.jet):
     '''
     plot the image2d
     
     :param vmin: minimum value for the colorbar
     :type vmin: float
     :param vmax: maximun value for the colorbar
     :type vmax: float
     :param colorbarcenter: do you want center the colorbar around 0
     :type colorbarcenter: bool
     :param colorbar: colorbar from matplotlib.cm
     
     .. note:: colorbar : cm.jet for eqstrain-stress
     '''
     if np.isnan(vmin):
         vmin=np.nanmin(self.field)
         
     if np.isnan(vmax):
         vmax=np.nanmax(self.field)
     
     # size of the image2d
     ss=np.shape(self.field)
     # create image
     img=plt.imshow(self.field,aspect='equal',extent=(0,ss[1]*self.res,0,ss[0]*self.res),cmap=colorbar,vmin=vmin,vmax=vmax)
     
     if colorbarcenter:
         zcolor=np.max(np.max(np.abs(self.field)))
         plt.clim(-zcolor, zcolor)
     
     # set up colorbar
     plt.colorbar(img,orientation='vertical',aspect=4)

    def plot_table_as_heatmap(self, intensity_min = 10**4, intensity_max = 10**9):
        table = self.get_table()
        # Numpy does not like dealing with NaN unless things are floats evidently.
        table_values = np.array(table.values, dtype=np.float64)
        # Since we will probably plot this with a log scale, let us set all NaN values to 1.
        table_values[np.isnan(table_values)] = 1

        # Do the plotting
        ski.io.imshow(table_values, norm=LogNorm(), cmap='cubehelix', interpolation='None')
        plt.clim(intensity_min, intensity_max)
        plt.grid(False)

        # Label the axis ticks so that they look like a plate

        ytick_positions = np.arange(0, self.well_rows)
        ytick_labels = ytick_positions + 1
        # We actually want the ytick labels to be letters
        ytick_labels = [chr(z + LOWERCASE_TO_INT).upper() for z in ytick_labels]

        plt.yticks(ytick_positions, ytick_labels)

        xtick_positions = np.arange(0, self.well_columns)
        xtick_labels = xtick_positions + 1

        plt.xticks(xtick_positions, xtick_labels)

def field(domain, z, title, clim = None,  path=None, save=True, show =
          False, ics=1, ext='.png', cmap=plt.cm.jet, fmt=None):
    """
    Given a domain, plot the nodal value z
   
    :param domain: :class:`~polyadcirc.run_framework.domain`
    :param z: :class:`numpy.ndarray`
    :param string title: plot title
    :param clim: :class:`numpy.clim`
    :type path: string or None
    :param path: directory to store plots
    :type save: bool
    :param save: flag for whether or not to save plots
    :type show: bool
    :param show: flag for whether or not to show plots
    :param int ics:  polar coordinate option (1 = cart coords, 2 = polar
        coords)
    :param string ext: file extension
    :param callable fmt: formatter for color bar, takes ``(x, pos)`` returns
        ``string`` 

    """
    if path is None:
        path = os.getcwd()
    plt.figure()
    plt.tripcolor(domain.triangulation, z, shading='gouraud',
                  cmap=cmap)
    plt.gca().set_aspect('equal')
    plt.autoscale(tight=True)
    if clim:
        plt.clim(clim[0], clim[1])
    add_2d_axes_labels(ics=ics)    
    colorbar(fmt=fmt)
    #plt.title(title)
    save_show(os.path.join(path, 'figs', title), save, show, ext)

def reduce_and_plot_clusters(X, model, title="", just_labelled=False, event_clusters=[]):

    if not type(model.labels_) is list:
        n_clus = len(set(model.labels_.tolist()))
        labels_ = model.labels_.tolist()
    else:
        labels_ = model.labels_
        n_clus = len(set(model.labels_))

    if just_labelled:
        idx = [i for i, label in enumerate(labels_) if label in event_clusters]
        print idx
        X = X[idx, :]
        labels_ = [labels_[i] for i in idx]
        label_dict = dict(zip(set(labels_), range(0, len(set(labels_)))))
        labels_ = [label_dict[label] for label in labels_]
        print labels_
        n_clus = len(event_clusters)
        print n_clus

    X_reduced = TruncatedSVD().fit_transform(X)
    X_embedded = TSNE(learning_rate=100).fit_transform(X_reduced)

    plt.figure(figsize=(10, 10))
    plt.title(title)
    plt.scatter(X_embedded[:, 0], X_embedded[:, 1], marker="x",
                c=labels_, cmap=plt.cm.get_cmap("jet", n_clus))
    plt.colorbar(ticks=min(range(500), range(n_clus)))
    plt.clim(-0.5, (n_clus - 0.5))
    plt.savefig('i/reduce_and_plot_clusters' +
                str(datetime.now()) + '.png')
    plt.show()

def plotColor(filename):
    fig = plt.figure(figsize=(10,20))
    plt.xlabel('$ x $')
    plt.ylabel('$ y $')
    z = F[0]
    im = plt.imshow(z, animated=True)
    plt.clim(0,0.5)
    plt.colorbar(im, label='wave function probability $ \parallel \psi ( x, y ) \parallel  ^ 2 $ ')


    def updatefig(i):
        try:
            im.set_array(F[i])
        except TypeError:
            quit()
        return im,

    # UNDER CONSTRUCTION
    # Set up formatting for the movie files
    #Writer = animation.writers['mencoder']
    #writer = Writer(fps=30, metadata=dict(artist='Me'), bitrate=1800)

    ani = animation.FuncAnimation(fig, updatefig, interval=speed, blit=True)
    # UNDER CONSTRUCTION
    #ani.save(filename + '.mp4', writer=writer)
    show()