def rootmov( numframes, degree, bins, dpi):
        #Main loop for making frame images
        for frame in range(1,numframes + 1):
               realy = list()
               imagy = list()
               percent = 1.0 * frame / numframes
               # Find the roots of all polynomials of given degree as coefficients vary.
	       for group in product(pathcoeff(percent),repeat=degree):
                        rootie = np.roots(group)
                        for rooter in list(rootie):
                                if rooter.imag != 0:
                                        realy.append(rooter.real)
                                        imagy.append(- rooter.imag)
               # Make histogram of roots.
               H, xedges, yedges = np.histogram2d(realy,imagy, bins=bins)
               H = np.log1p( 1 / (1 + H ) )
               # Configure and save an image of the histogram.
               fig=plt.figure( facecolor='k', edgecolor='k')
               ax=plt.gca()
               plt.setp(ax, frame_on=True)
               plt.setp(ax.get_xticklabels(), visible=False)
               plt.setp(ax.get_yticklabels(), visible=False)
               plt.setp(ax.get_xticklines(), visible=False)
               plt.setp(ax.get_yticklines(), visible=False)
               plt.imshow(H,interpolation='bicubic',extent=[0,1000,0,600], cmap=dynacm( percent ) )
               plt.savefig("root_test{:04}.png".format(frame),dpi=dpi, facecolor='k', edgecolor='k', bbox_inches='tight')
               ax.clear()
               plt.close(fig)

def main():
    gw = gridworld()
    a = agent(gw)

    for epoch in range(20):
        a.initEpoch()
        while True:
            rwd, stat, act = a.takeAction()
            a.updateQ(rwd, stat, act)
            if gw.status() == 'Goal':
                break
            if mod(a.counter, 10)==0:
                print(gw.state())
                print(gw.field())
        print('Finished')
        print(a.counter)
        print(gw.state())
        print(gw.field())
        Q = transpose(a.Q(), (2,0,1))
        for i in range(4):
            plt.subplot(2,2,i)
            plt.imshow(Q[i], interpolation='nearest')
            plt.title(a.actions()[i])
            plt.colorbar()
        plt.show()

def show_overlay(img3d, cc3d, ncc=10, s=85, xyz = 'xy',alpha=.8):
    """Shows the connected components overlayed over img3d

    Input
    ======
    img3d -- 3d array
    cc3d -- 3d array ( preferably of same shape as img3d, use get_3d_cc(...) )
    ncc -- where to cut off the color scale
    s -- slice to show
    xyz -- which projection to use in {'xy','xz','yz'}
    """
    cc = get_slice(cc3d,s,xyz)
    img = get_slice(img3d,s,xyz)

    notcc = np.isnan(cc)
    incc = np.not_equal(notcc,True)

    img4 = plt.cm.gray(img/np.nanmax(img))
    if ncc is not np.Inf:
        cc = plt.cm.jet(cc/float(ncc))
    else:
        cc = plt.cm.jet(np.log(cc)/np.log(np.nanmax(cc)))

    cc[notcc,:]=img4[notcc,:]
    cc[incc,3] = 1-img[incc]/(2*np.nanmax(img))

    plt.imshow(cc)

def export(data, F, k):
    '''Write data to a png image
    
    Arguments
    ---------
    data : numpy.ndarray
        array containing the data to be written as png image
    F : float
        feed rate of the current configuration
    k : float
        rate constant of the current configuration
    '''
        
    figsize = tuple(s / 72.0 for s in data.shape)
    fig = plt.figure(figsize=figsize, dpi=72.0, facecolor='white')
    fig.add_axes([0, 0, 1, 1], frameon=False)
    plt.xticks([])
    plt.yticks([])

    plt.imshow(data, cmap=plt.cm.RdBu_r, interpolation='bicubic')
    plt.gci().set_clim(0, 1)

    filename = './study/F{:03d}-k{:03d}.png'.format(int(1000*F), int(1000*k))
    plt.savefig(filename, dpi=72.0)
    plt.close()

    def Test(self):
        test_Dir = "Result";        
        if not os.path.exists(test_Dir):
            os.makedirs(test_Dir);

        test_Label_List = [0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5];
        test_Label_Pattern = np.zeros((16, 10));
        test_Label_Pattern[np.arange(16), test_Label_List] = 1;            
        feed_Dict = {
            self.noise_Placeholder: np.random.uniform(-1., 1., size=[16, self.noise_Size]),
            self.label_for_Fake_Placeholder: test_Label_Pattern,
            self.is_Training_Placeholder: False
            };   #Batch is constant in the test.
        global_Step, mnist_List = self.tf_Session.run(self.test_Tensor_List, feed_dict = feed_Dict);

        fig = plt.figure(figsize=(4, 4))
        gs = gridspec.GridSpec(4, 4)
        gs.update(wspace=0.05, hspace=0.05)

        for index, mnist in enumerate(mnist_List):
            ax = plt.subplot(gs[index])
            plt.axis('off')
            ax.set_xticklabels([])
            ax.set_yticklabels([])
            ax.set_aspect('equal')
            plt.imshow(mnist.reshape(28, 28), cmap='Greys_r')

        plt.savefig('%s/S%d.png' % (test_Dir, global_Step), bbox_inches='tight');
        plt.close();

    def atest_interpolation_coast(self):
        """Test interpolation."""
        reader = reader_ROMS_native.Reader('/disk2/data/SVIM/ocean_avg_20081201.nc')
        num_points = 50
        np.random.seed(0)  # To get the same random numbers each time
        lons = np.random.uniform(12, 16, num_points)
        lats = np.random.uniform(68.3, 68.3, num_points)
        z = np.random.uniform(-100, 0, num_points)
        x, y = reader.lonlat2xy(lons, lats)

        variables = ['x_sea_water_velocity', 'y_sea_water_velocity',
                     'sea_water_temperature']
        # Read a block of data covering the points
        data = reader.get_variables(variables, time=reader.start_time,
                                    x=x, y=y, z=z, block=True)
        import matplotlib.pyplot as plt
        plt.imshow(data['x_sea_water_velocity'][0,:,:])
        plt.colorbar()
        plt.show()

        b = ReaderBlock(data, interpolation_horizontal='nearest')

        env, prof = b.interpolate(x, y, z, variables,
                                  profiles=['x_sea_water_velocity'],
                                  profiles_depth=[-30, 0])

def heatmap(vals, size=6, aspect=1):
    """
    Plot a heatmap from matrix data
    """
    plt.figure(figsize=(size, size))
    plt.imshow(vals, cmap="gray", aspect=aspect, interpolation="none", vmin=0, vmax=1)
    plt.axis("off")

def plot_images(data_list, data_shape="auto", fig_shape="auto"):
    """
    plotting data on current plt object.
    In default,data_shape and fig_shape are auto.
    It means considered the data as a sqare structure.
    """
    n_data = len(data_list)
    if data_shape == "auto":
        sqr = int(n_data ** 0.5)
        if sqr * sqr != n_data:
            data_shape = (sqr + 1, sqr + 1)
        else:
            data_shape = (sqr, sqr)
    plt.figure(figsize=data_shape)

    for i, data in enumerate(data_list):
        plt.subplot(data_shape[0], data_shape[1], i + 1)
        plt.gray()
        if fig_shape == "auto":
            fig_size = int(len(data) ** 0.5)
            if fig_size ** 2 != len(data):
                fig_shape = (fig_size + 1, fig_size + 1)
            else:
                fig_shape = (fig_size, fig_size)
        Z = data.reshape(fig_shape[0], fig_shape[1])
        plt.imshow(Z, interpolation="nearest")
        plt.tick_params(labelleft="off", labelbottom="off")
        plt.tick_params(axis="both", which="both", left="off", bottom="off", right="off", top="off")
        plt.subplots_adjust(hspace=0.05)
        plt.subplots_adjust(wspace=0.05)

def draw_digit(data):
    size = int(len(data) ** 0.5)
    Z = data.reshape(size, size)  # convert from vector to matrix
    plt.imshow(Z, interpolation="None")
    plt.gray()
    plt.tick_params(labelbottom="off")
    plt.tick_params(labelleft="off")

    def psf_dl(self, 
        N_OS = None,
        plate_scale_as_px = None,
        plotit = False,
        return_efield = False
        ):
        """  Find the PSF of the wavefront at a given wavelength. 

        Parameters
        ----------
        return_efield: boolean
            Do we return an electric field? If not, return the intensity. """

        # Make the field uniform.
        self.flatten_field()

        # The PSF is then simply the FFT of the pupil.
        psf = self.image(N_OS = N_OS, plate_scale_as_px = plate_scale_as_px, 
            return_efield = return_efield)
        psf /= sum(psf.flatten())
        if plotit:
            axesScale = [0, self.sz*self.m_per_px, 0, self.sz*self.m_per_px]
            plt.figure()
            plt.imshow(psf, extent = axesScale)
            plt.title('PSF of optical system')
        return psf
            

            

def template_matching():
    img = cv2.imread('messi.jpg',0)
    img2 = img.copy()
    template = cv2.imread('face.png',0)
    w, h = template.shape[::-1]

    # All the 6 methods for comparison in a list
    methods = ['cv2.TM_CCOEFF', 'cv2.TM_CCOEFF_NORMED', 'cv2.TM_CCORR',
            'cv2.TM_CCORR_NORMED', 'cv2.TM_SQDIFF', 'cv2.TM_SQDIFF_NORMED']

    for meth in methods:
        img = img2.copy()
        method = eval(meth)

        # Apply template Matching
        res = cv2.matchTemplate(img,template,method)
        min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(res)

        # If the method is TM_SQDIFF or TM_SQDIFF_NORMED, take minimum
        if method in [cv2.TM_SQDIFF, cv2.TM_SQDIFF_NORMED]:
            top_left = min_loc
        else:
            top_left = max_loc
        bottom_right = (top_left[0] + w, top_left[1] + h)

        cv2.rectangle(img,top_left, bottom_right, 255, 2)

        plt.subplot(121),plt.imshow(res,cmap = 'gray')
        plt.title('Matching Result'), plt.xticks([]), plt.yticks([])
        plt.subplot(122),plt.imshow(img,cmap = 'gray')
        plt.title('Detected Point'), plt.xticks([]), plt.yticks([])
        plt.suptitle(meth)

        plt.show()

def plot_confusion_matrix(cm, classes,
                          normalize=False,
                          title='Confusion matrix',
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=45)
    plt.yticks(tick_marks, classes)

    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
    else:
        print('Confusion matrix, without normalization')

    print(cm)

    thresh = cm.max() / 2.
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        plt.text(j, i, cm[i, j],
                 horizontalalignment="center",
                 color="white" if cm[i, j] > thresh else "black")

    plt.tight_layout()
    plt.ylabel('True label')
    plt.xlabel('Predicted label')

def subaps_to_grid(subap_frame, plot=False):
    """
    Turn a 1D list of subaperture values to an 11x11 grid
    
    See figure from Rudy / Contreras for how we count our subaps.
    
    Or just do this:
    
    >>> subaps_to_grid(range(97), plot=True)
    >>> overlay_indices()
    
    (n.b. sensible image plotting requires the
    origin='bottom' keyword argument unless you set
    it as default in matplotlibrc)
    """
    grid = np.ndarray((11,11))
    grid[:,:] = np.nan
    grid[0][3:8] = subap_frame[0:5]
    grid[1][2:9] = subap_frame[5:12]
    grid[2][1:10] = subap_frame[12:21]
    grid[3] = subap_frame[21:32]
    grid[4] = subap_frame[32:43]
    grid[5] = subap_frame[43:54]
    grid[6] = subap_frame[54:65]
    grid[7] = subap_frame[65:76]
    grid[8][1:10] = subap_frame[76:85]
    grid[9][2:9] = subap_frame[85:92]
    grid[10][3:8] = subap_frame[92:97]
    grid = grid.T # we're filling in row-by-row from the top, but numbering starts
                  # in the bottom left with zero and proceeds column-by-column
    if plot:
        plt.imshow(grid, origin='bottom')
    return grid

def imview(*args, **kwargs):
    """
    A more sensible matplotlib-based image viewer command,
    a wrapper around `matplotlib.pyplot.imshow`.

    Parameters are identical to `matplotlib.pyplot.imshow`
    but this behaves somewhat differently:

      * By default, it creates a new figure (unless a
        `figure` keyword argument is supplied.
      * It modifies the axes of that figure to use the
        full frame, without ticks or tick labels.
      * It turns on `nearest` interpolation by default
        (i.e., it does not antialias pixel data). This
        can be overridden with the `interpolation`
        argument as in `imshow`.

    All other arguments and keyword arguments are passed
    on to `imshow`.`
    """
    if 'figure' not in kwargs:
        f = plt.figure()
    else:
        f = kwargs['figure']
    new_ax = matplotlib.axes.Axes(f, [0, 0, 1, 1],
                                  xticks=[], yticks=[],
                                  frame_on=False)
    f.delaxes(f.gca())
    f.add_axes(new_ax)
    if len(args) < 5 and 'interpolation' not in kwargs:
        kwargs['interpolation'] = 'nearest'
    plt.imshow(*args, **kwargs)

def plot_bold_nii(filename, timepoint):
    """Plot all slices of a fMRI image in one plot at a specific time point

    Parameters:
    -----------
    filename: BOLD.nii.gz
    timepoint: the time point chose

    Return:
    -------
    None

    Note:
    -----
    The function produce a plot
    """
    img = nib.load(filename)
    data = img.get_data()
    assert timepoint <= data.shape[-1]
    plot_per_row = int(np.ceil(np.sqrt(data.shape[2])))
    frame = np.zeros((data.shape[0]*plot_per_row, data.shape[1]*plot_per_row))
    num_of_plots = 0
    for i in range(plot_per_row):
        j = 0
        while j < plot_per_row and num_of_plots < data.shape[2]:
            frame[data.shape[0]*i:data.shape[0]*(i+1), data.shape[1]*j:data.shape[1]*(j+1)] = data[:,:,num_of_plots,timepoint]
            num_of_plots+=1
            j+=1
    plt.imshow(frame, interpolation="nearest",cmap="gray")
    return None

	def logLikelihood(self, data):
		assert self.numImages == data.numImages

		which = [np.nonzero(data.id == i)[0]\
				for i in xrange(0, self.numImages)]

		# Mean vector
		m = np.empty(data.t.size)
		for i in xrange(0, self.numImages):
			m[which[i]] = self.mag[i]

		# Covariance matrix
		[t1, t2] = np.meshgrid(data.t, data.t)
		lags = np.abs(t2 - t1)
		plt.imshow(lags)
		plt.show()
		ids = np.meshgrid(data.id, data.id)


		equal = ids[0] == ids[1]

##		try:
##			L = la.cholesky(C)
##		except:
##			return -np.inf
#		y = data.y - m
#		logDeterminant = 2.0*np.sum(np.log(np.diag(L)))
#		solution = la.cho_solve((L, True), y)
#		exponent = np.dot(y, solution)
#		logL = -0.5*data.t.size*np.log(2.0*np.pi) - 0.5*logDeterminant - 0.5*exponent
		return 0.

def plot_jacobian(A, name, cmap= plt.cm.coolwarm, normalize=True, precision=1e-6):

    """
    Customized visualization of jacobian matrices for observing
    sparsity patterns
    """
    
    plt.figure()
    fig, ax = plt.subplots()
    
    if normalize is True:
        plt.imshow(A, interpolation='none', cmap=cmap,
                   norm = mpl.colors.Normalize(vmin=-1.,vmax=1.))
    else:
        plt.imshow(A, interpolation='none', cmap=cmap)        
    plt.colorbar(format=ticker.FuncFormatter(fmt))
    
    ax.spy(A, marker='.', markersize=0,  precision=precision)
    
    ax.spines['right'].set_visible(True)
    ax.spines['bottom'].set_visible(True)
    ax.xaxis.set_ticks_position('top')
    ax.yaxis.set_ticks_position('left')

    xlabels = np.linspace(0, A.shape[0], 5, True, dtype=int)
    ylabels = np.linspace(0, A.shape[1], 5, True, dtype=int)

    plt.xticks(xlabels)
    plt.yticks(ylabels)

    plt.savefig(name, bbox_inches='tight', pad_inches=0.05)
    
    plt.close()

    return

def draw():
    global g_matrix
    m = np.matrix(np.pad(g_matrix,((1,1),(1,1)), mode='constant', constant_values=0))
    plt.imshow(m, cmap = cm.Greys_r, interpolation='nearest')
    # 横纵坐标
    # plt.xticks([]), plt.yticks([])
    plt.show()

def edge_detect(img):
    BLUR_SIZE = 51
    TRUNC_RATIO = 0.75
    CLOSING_SIZE = 5

    # denoised = cv2.fastNlMeansDenoisingColored(img,None,10,10,7,21)
    # img = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    # too_bright=np.logical_and(img[:,:,1]<50, img[:,:,2]>200)
    # np.set_printoptions(threshold=np.nan)
    # np.savetxt('conconcon',img[:,:,1],'%i')
    # img[:,:,1]=np.where(too_bright, np.sqrt(img[:,:,1])+70, img[:,:,1])
    # img = cv2.cvtColor(img, cv2.COLOR_HSV2BGR)

    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    blur = cv2.blur(gray, (BLUR_SIZE, BLUR_SIZE))

    edge = np.floor(0.5 * gray + 0.5 * (255 - blur)).astype('uint8')

    hist,bins = np.histogram(edge.flatten(), 256, [0, 256])
    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max() / cdf.max()
    cdf_m = np.ma.masked_equal(cdf, 0)
    cdf_m = (cdf_m - cdf_m.min()) * 255 / (cdf_m.max() - cdf_m.min())
    cdf = np.ma.filled(cdf_m, 0).astype('uint8')
    equ = cdf[edge]

    hist,bins = np.histogram(equ.flatten(),256,[0,256])
    max_idx = np.argmax(hist);
    hist_clean = np.where(equ > TRUNC_RATIO * max_idx, 255, equ)

    kernel = np.ones((CLOSING_SIZE, CLOSING_SIZE), np.uint8)
    closing = cv2.morphologyEx(hist_clean, cv2.MORPH_CLOSE, kernel)
    plt.imshow(closing, cmap='Greys_r')
    plt.show()
    cv2.waitKey(100)

def plot_brights(ax, path, star, regionList, goal=False):
    '''
    Components of this routine:
        Projected brightness map
         
    Please note that this has been modified for use in diagnostic plots, 
    there should really be a way to specify a windowNumber for real data
    '''
    currentWindow = 0

    ###########################
    # Make the brightness map #
    ###########################
    img = make_bright_image(star, regionList, currentWindow, goal=goal)
    
    plt.imsave(path + "temp.jpg", img, cmap='hot', vmin=0.85, vmax=1.15)
    plt.imshow(img, cmap='hot')
    #Create the plot
    bmap = Basemap(projection='moll', lon_0 = 0, ax=ax)
    bmap.warpimage(path + "temp.jpg", ax=ax)
    
    if goal:
        ax.set_title("Desired Map")
    else:
        ax.set_title("Average Map")