def plot_tracks(src, fakewcs, spa=None, **kwargs):
    # NOTE -- MAGIC 61 = monthly; this is ASSUMEd below.
    tt = np.linspace(2010., 2015., 61)
    t0 = TAITime(None, mjd=TAITime.mjd2k + 365.25*10)
    #rd0 = src.getPositionAtTime(t0)
    #print 'rd0:', rd0
    xx,yy = [],[]
    rr,dd = [],[]
    for t in tt:
        #print 'Time', t
        rd = src.getPositionAtTime(t0 + (t - 2010.)*365.25*24.*3600.)
        ra,dec = rd.ra, rd.dec
        rr.append(ra)
        dd.append(dec)
        ok,x,y = fakewcs.radec2pixelxy(ra,dec)
        xx.append(x - 1.)
        yy.append(y - 1.)

    if spa is None:
        spa = [None,None,None]
    for rows,cols,sub in spa:
        if sub is not None:
            plt.subplot(rows,cols,sub)
        ax = plt.axis()
        plt.plot(xx, yy, 'k-', **kwargs)
        plt.axis(ax)

    return rr,dd,tt

def main():
    SAMPLE_NUM = 10
    degree = 9
    x, y = sin_wgn_sample(SAMPLE_NUM)
    fig = pylab.figure(1)
    pylab.grid(True)
    pylab.xlabel('x')
    pylab.ylabel('y')
    pylab.axis([-0.1,1.1,-1.5,1.5])

    # sin(x) + noise
    # markeredgewidth mew
    # markeredgecolor mec
    # markerfacecolor mfc

    # markersize      ms
    # linewidth       lw
    # linestyle       ls
    pylab.plot(x, y,'bo',mew=2,mec='b',mfc='none',ms=8)

    # sin(x)
    x2 = linspace(0, 1, 1000)
    pylab.plot(x2,sin(2*x2*pi),'#00FF00',lw=2,label='$y = \sin(x)$')

    # polynomial fit
    reg = exp(-18)
    w = curve_poly_fit(x, y, degree,reg) #w = polyfit(x, y, 3)
    po = poly1d(w)      
    xx = linspace(0, 1, 1000)
    pylab.plot(xx, po(xx),'-r',label='$M = 9, \ln\lambda = -18$',lw=2)
    
    pylab.legend()
    pylab.show()
    fig.savefig("poly_fit9_10_reg.pdf")

def plotLDDecaySelection3d(ax, sweep=False):
    import pylab as plt; import matplotlib as mpl;mpl.rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size':16}) ;    mpl.rc('text', usetex=True)

    def neutral(ld0, t, d, r=2 * 1e-8):
        if abs(d) <= 5e3:
            d = np.sign(d) * 5e3
        if d == 0:
            d = 5e3
        return ((np.exp(-2 * r * t * abs(d)))) * ld0

    t = np.arange(0, 200 + 1., 2)
    L=1e6+1
    pos=500000
    r=2*1e-8
    ld0 = 0.5
    s = 0.05
    nu0 = 0.1
    positions=np.arange(0,L,1000)
    dist=(positions - pos)
    T, D = np.meshgrid(t, dist)
    if not sweep:
        zs = np.array([neutral(ld0, t, d) for t, d in zip(np.ravel(T), np.ravel(D))])
    else:
        zs = np.array([LD(t, ld0, s, nu0, r, abs(d), 0) for t, d in zip(np.ravel(T), np.ravel(D))])
    Z = zs.reshape(T.shape)
    ax.plot_surface(T, D, Z,cmap=mpl.cm.autumn)
    ax.set_xlabel('Generations')
    ax.set_ylabel('Position')
    plt.yticks(plt.yticks()[0][1:-1],map(lambda x:'{:.0f}K'.format((pos+(x))/1000),plt.yticks()[0][1:-1]))
    plt.ylim([-500000,500000])
    ax.set_zlabel(r"$|\rho_t|$")
    pplt.setSize(plt.gca(), fontsize=6)
    plt.axis('tight');

def compareAnimals(animals, precision):
    """Assumes animals is a list of animals, precision an int >= 0
       Builds a table of Euclidean distance between each animal"""
    #Get labels for columns and rows
    columnLabels = []
    for a in animals:
        columnLabels.append(a.getName())
    rowLabels = columnLabels[:]
    tableVals = []
    #Get distances between pairs of animals
    #For each row
    for a1 in animals:
        row = []
        #For each column
        for a2 in animals:
            if a1 == a2:
                row.append('--')
            else:
                distance = a1.distance(a2)
                row.append(str(round(distance, precision)))
        tableVals.append(row)
    #Produce table
    table = pylab.table(rowLabels = rowLabels,
                        colLabels = columnLabels,
                        cellText = tableVals,
                        cellLoc = 'center',
                        loc = 'center',
                        colWidths = [0.2]*len(animals))
    table.scale(1, 2.5)
    pylab.axis('off') #Don't display x and y-axes
    pylab.savefig('distances')

def display_image_from_array(nparray,colory='binary',roi=None):
	"""
	Produce a display of the nparray 2D matrix
	@param nparray : image to display
	@type nparray : numpy 2darray
	@param colory : color mapping of the image (see http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps)
	@type colory : string
	"""
	#Set the region of interest to display :
	#  (0,0) is set at lower left corner of the image
	if roi == None:
		roi = ((0,0),nparray.shape)
		nparraydsp = nparray
		print roi
	elif type(roi[0])==tuple and type(roi[1])==tuple: 
		# Case of 2 points definition of the domain : roi = integers index of points ((x1,y1),(x2,y2))
		print roi
		nparraydsp = nparray[roi[0][0]:roi[1][0],roi[0][1]:roi[1][1]]
	elif type(roi[0])==int and type(roi[1])==int:	
		# Case of image centered domain : roi = integers (width,high)
		nparraydsp = nparray[int(nparray.shape[0]/2)-int(roi[0])/2:int(nparray.shape[0]/2)+int(roi[0])/2,int(nparray.shape[1]/2)-int(roi[1])/2:int(nparray.shape[1]/2)+int(roi[1])/2]
	fig = pylab.figure()
    #Display array with grayscale intensity and no pixel smoothing interpolation
	pylab.imshow(nparraydsp,cmap=colory,interpolation='nearest')#,origin='lower')
	pylab.colorbar()
	pylab.axis('off')

	def update(self):
		if self.pose != []:
			plt.figure(1)
			clf()
			self.fig1 = plt.figure(num=1, figsize=(self.window_size, \
				self.window_size), dpi=80, facecolor='w', edgecolor='w')
			title (self.title)			
			xlabel('Easting [m]')
			ylabel('Northing [m]')
			axis('equal')
			grid (True)
			poseT = zip(*self.pose)	
			pose_plt = plot(poseT[1],poseT[2],'#ff0000')

			if self.wptnav != []:
				mode = self.wptnav[-1][MODE]

				if not (self.wptnav[-1][B_E] == 0 and self.wptnav[-1][B_N] == 0 and self.wptnav[-1][A_E] == 0 and self.wptnav[-1][A_N] == 0):
					b_dot = plot(self.wptnav[-1][B_E],self.wptnav[-1][B_N],'ro',markersize=8)
					a_dot = plot(self.wptnav[-1][A_E],self.wptnav[-1][A_N],'go',markersize=8)
					ab_line = plot([self.wptnav[-1][B_E],self.wptnav[-1][A_E]],[self.wptnav[-1][B_N],self.wptnav[-1][A_N]],'g')
					target_dot = plot(self.wptnav[-1][TARGET_E],self.wptnav[-1][TARGET_N],'ro',markersize=5)

				if mode == -1:
					pose_dot = plot(self.wptnav[-1][POSE_E],self.wptnav[-1][POSE_N],'b^',markersize=8)
				elif mode == 1:
					pose_dot = plot(self.wptnav[-1][POSE_E],self.wptnav[-1][POSE_N],'bs',markersize=8)
				elif mode == 2:
					pose_dot = plot(self.wptnav[-1][POSE_E],self.wptnav[-1][POSE_N],'bo',markersize=8)

		if self.save_images:
			self.fig1.savefig ('plot_map%05d.jpg' % self.image_count)
			self.image_count += 1
		draw()

def plotslice(pos,filename='',boxsize=100.):
    ng = pos.shape[0]
    M.clf()
    M.scatter(pos[ng/4,:,:,1].flatten(),pos[ng/4,:,:,2].flatten(),s=1.,lw=0.)
    M.axis('tight')
    if filename != '':
        M.savefig(filename)

def plot_iris_knn():
    iris = datasets.load_iris()
    X = iris.data[:, :2]  # we only take the first two features. We could
                        # avoid this ugly slicing by using a two-dim dataset
    y = iris.target

    knn = neighbors.KNeighborsClassifier(n_neighbors=3)
    knn.fit(X, y)

    x_min, x_max = X[:, 0].min() - .1, X[:, 0].max() + .1
    y_min, y_max = X[:, 1].min() - .1, X[:, 1].max() + .1
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),
                         np.linspace(y_min, y_max, 100))
    Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pl.figure()
    pl.pcolormesh(xx, yy, Z, cmap=cmap_light)

    # Plot also the training points
    pl.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold)
    pl.xlabel('sepal length (cm)')
    pl.ylabel('sepal width (cm)')
    pl.axis('tight')

def plot_polynomial_regression():
    rng = np.random.RandomState(0)
    x = 2*rng.rand(100) - 1
    f = lambda t: 1.2 * t**2 + .1 * t**3 - .4 * t **5 - .5 * t ** 9
    y = f(x) + .4 * rng.normal(size=100)

    x_test = np.linspace(-1, 1, 100)

    pl.figure()
    pl.scatter(x, y, s=4)

    X = np.array([x**i for i in range(5)]).T
    X_test = np.array([x_test**i for i in range(5)]).T
    regr = linear_model.LinearRegression()
    regr.fit(X, y)
    pl.plot(x_test, regr.predict(X_test), label='4th order')

    X = np.array([x**i for i in range(10)]).T
    X_test = np.array([x_test**i for i in range(10)]).T
    regr = linear_model.LinearRegression()
    regr.fit(X, y)
    pl.plot(x_test, regr.predict(X_test), label='9th order')

    pl.legend(loc='best')
    pl.axis('tight')
    pl.title('Fitting a 4th and a 9th order polynomial')

    pl.figure()
    pl.scatter(x, y, s=4)
    pl.plot(x_test, f(x_test), label="truth")
    pl.axis('tight')
    pl.title('Ground truth (9th order polynomial)')

def plot_adsorbed_circles(adsorbed_x,adsorbed_y,radius, width, label=""):
    import pylab
    from matplotlib.patches import Circle

    # Plot each run
    fig = pylab.figure()
    ax = fig.add_subplot(111)
    for p in range(len(adsorbed_x)):
        ax.add_patch(Circle((adsorbed_x[p], adsorbed_y[p]), radius))
        # Plot "image" particles to verify that periodic boundary conditions are working
#        if adsorbed_x[p] < radius:
#            ax.add_patch(Circle((adsorbed_x[p] + width,adsorbed_y[p]), radius, facecolor='red'))
#        elif adsorbed_x[p] > (width-radius):
#            ax.add_patch(Circle((adsorbed_x[p] - width,adsorbed_y[p]), radius, facecolor='red'))
#        if adsorbed_y[p] < radius:
#            ax.add_patch(Circle((adsorbed_x[p],adsorbed_y[p] + width), radius, facecolor='red'))
#        elif adsorbed_y[p] > (width-radius):
#            ax.add_patch(Circle((adsorbed_x[p],adsorbed_y[p] - width), radius, facecolor='red'))

    ax.set_aspect(1.0)
    pylab.axhline(y=0, color='k')
    pylab.axhline(y=width, color='k')
    pylab.axvline(x=0, color='k')
    pylab.axvline(x=width, color='k')
    pylab.axis([-0.1*width, width*1.1, -0.1*width, width*1.1])
    pylab.xlabel("non-dimensional x")
    pylab.ylabel("non-dimensional y")
    pylab.title("Adsorbed particles at theta="+label)

    return ax

def hinton(W, maxWeight=None):
    """
    Source: http://wiki.scipy.org/Cookbook/Matplotlib/HintonDiagrams
    Draws a Hinton diagram for visualizing a weight matrix.
    Temporarily disables matplotlib interactive mode if it is on,
    otherwise this takes forever.
    """
    reenable = False
    if pl.isinteractive():
        pl.ioff()
    pl.clf()
    height, width = W.shape
    if not maxWeight:
        maxWeight = 2**np.ceil(np.log(np.max(np.abs(W)))/np.log(2))

    pl.fill(np.array([0,width,width,0]),np.array([0,0,height,height]),'gray')
    pl.axis('off')
    pl.axis('equal')
    for x in xrange(width):
        for y in xrange(height):
            _x = x+1
            _y = y+1
            w = W[y,x]
            if w > 0:
                _blob(_x - 0.5, height - _y + 0.5, min(1,w/maxWeight),'white')
            elif w < 0:
                _blob(_x - 0.5, height - _y + 0.5, min(1,-w/maxWeight),'black')
    if reenable:
        pl.ion()
    pl.show()

def draw_io( type ):
    f = open("data\\io.txt")
    s_list = f.readlines()
    f.close()

    r_list_x = []
    r_list_y = []
    w_list_x = []
    w_list_y = []

    y_min = 0x80000000
    y_max = 0

    for s in s_list:
        pos = s.find('[')
        addr = int( s[pos+1:pos+9], 16 )
        if addr > y_max: y_max = addr
        if addr < y_min: y_min = addr 
        if s.find('|R') != -1:
            r_list_y.append( addr )
            r_list_x.append( int( s.strip('\n').split('#')[1] ) ) #read counter
        if s.find('|W') != -1:
            w_list_y.append( int( addr ) )
            w_list_x.append( int( s.strip('\n').split('#')[1] ) ) #read counter

    if type == 'W': pylab.plot( w_list_x, w_list_y, "ro" )
    if type == 'R': pylab.plot( r_list_x, r_list_y, "ro" )

    pylab.axis( [0, get_trace_count(), y_min - 1000, y_max + 1000] )
    pylab.show()

    def __init__(self, baseConfig):
        self.figsize = baseConfig.get("figsize", None)
        self.axis = baseConfig.get("axis", None)
        self.title = baseConfig.get("title", "NoName")
        self.ylabel = baseConfig.get("ylabel", "NoName")
        self.grid = baseConfig.get("grid", False)
        self.xaxis_locator = baseConfig.get("xaxis_locator", None)
        self.yaxis_locator = baseConfig.get("yaxis_locator", None)
        self.legend_loc = baseConfig.get("legend_loc", 0)

        if self.figsize != None:
            pylab.figure(figsize=self.figsize)
        if self.axis != None:
            pylab.axis(self.axis)

        pylab.title(self.title)
        pylab.ylabel(self.ylabel)
        ax = pylab.gca()
        pylab.grid(self.grid)
        if self.xaxis_locator != None:
            ax.xaxis.set_major_locator(pylab.MultipleLocator(self.xaxis_locator))
        if self.yaxis_locator != None:
            ax.yaxis.set_major_locator(pylab.MultipleLocator(self.yaxis_locator))
        self.lineList = []
        self.id = 1

 def __init__(self, baseConfig) :
     self.figsize = baseConfig.get('figsize',None)
     self.axis = baseConfig.get('axis',None)
     self.title = baseConfig.get('title','NoName')
     self.ylabel = baseConfig.get('ylabel','NoName')
     self.grid = baseConfig.get('grid',False)
     self.xaxis_locator = baseConfig.get('xaxis_locator',None)
     self.yaxis_locator = baseConfig.get('yaxis_locator',None)
     self.legend_loc = baseConfig.get('legend_loc',0)
     
     if self.figsize != None :
         pylab.figure(figsize = self.figsize)
     if self.axis != None :
         pylab.axis(self.axis)
     
     pylab.title(self.title)
     pylab.ylabel(self.ylabel)
     ax = pylab.gca()
     pylab.grid(self.grid)
     if self.xaxis_locator != None :
         ax.xaxis.set_major_locator( pylab.MultipleLocator(self.xaxis_locator) )
     if self.yaxis_locator != None :
         ax.yaxis.set_major_locator( pylab.MultipleLocator(self.yaxis_locator) )
     self.lineList = []
     self.id = 1

    def plot_density(self, plot_filename="out/density.png"):
        x, y, labels = self.load_data()

        figure(figsize=(self.fig_width, self.fig_height), dpi=80)
        # Perform a kernel density estimator on the coords in data.
        # The following 10 lines can be commented out if density map not needed.
        space_factor = 1.2
        xmin = space_factor * x.min()
        xmax = space_factor * x.max()
        ymin = space_factor * y.min()
        ymax = space_factor * y.max()
        X, Y = mgrid[xmin:xmax:100j, ymin:ymax:100j]
        positions = c_[X.ravel(), Y.ravel()]
        values = c_[x, y]
        kernel = stats.kde.gaussian_kde(values.T)
        Z = reshape(kernel(positions.T).T, X.T.shape)
        imshow(rot90(Z), cmap=cm.gist_earth_r, extent=[xmin, xmax, ymin, ymax])

        # Plot the labels
        num_labels_to_plot = min([len(labels), self.max_labels, len(x), len(y)])
        if self.has_labels:
            for i in range(num_labels_to_plot):
                text(x[i], y[i], labels[i])  # assumes m size and order matches labels
        else:
            plot(x, y, "k.", markersize=1)
        axis("equal")
        axis("off")
        savefig(plot_filename)
        print "wrote %s" % (plot_filename)

    def plot_margin(X1_train, X2_train, clf):
        def f(x, w, b, c=0):
            # given x, return y such that [x,y] in on the line
            # w.x + b = c
            return (-w[0] * x - b + c) / w[1]

        pl.plot(X1_train[:,0], X1_train[:,1], "ro")
        pl.plot(X2_train[:,0], X2_train[:,1], "bo")
        pl.scatter(clf.sv[:,0], clf.sv[:,1], s=100, c="g")

        # w.x + b = 0
        a0 = -4; a1 = f(a0, clf.w, clf.b)
        b0 = 4; b1 = f(b0, clf.w, clf.b)
        pl.plot([a0,b0], [a1,b1], "k")

        # w.x + b = 1
        a0 = -4; a1 = f(a0, clf.w, clf.b, 1)
        b0 = 4; b1 = f(b0, clf.w, clf.b, 1)
        pl.plot([a0,b0], [a1,b1], "k--")

        # w.x + b = -1
        a0 = -4; a1 = f(a0, clf.w, clf.b, -1)
        b0 = 4; b1 = f(b0, clf.w, clf.b, -1)
        pl.plot([a0,b0], [a1,b1], "k--")

        pl.axis("tight")
        pl.show()

    def __call__(self, n):
        if len(self.f.shape) == 3:
            # f = f[x,v,t], 2 dim in phase space
            ft = self.f[n,:,:]
            pylab.pcolormesh(self.X, self.V, ft.T, cmap = 'jet')
            pylab.colorbar()
            pylab.clim(0,0.38) # for Landau test case
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$v$', fontsize = 18)
            pylab.title('$N_x$ = %d, $N_v$ = %d, $t$ = %2.1f' % (self.x.N, self.v.N, self.it*self.t.width))
            pylab.savefig(self.path + self.filename)
            pylab.clf()
            return None

        if len(self.f.shape) == 2:
            # f = f[x], 1 dim in phase space
            ft = self.f[n,:]
            pylab.plot(self.x.gridvalues,ft,'ob')
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$f(x)$', fontsize = 18)
            pylab.savefig(self.path + self.filename)
            return None

def show_frequencies(video_filename, bounds=None):
    """Graph the average value of the video as well as the frequency strength"""
    original_video, fps = load_video(video_filename)
    print fps
    averages = []

    if bounds:
        for x in range(1, original_video.shape[0] - 1):
            averages.append(original_video[x, bounds[2]:bounds[3], bounds[0]:bounds[1], :].sum())
    else:
        for x in range(1, original_video.shape[0] - 1):
            averages.append(original_video[x, :, :, :].sum())

    charts_x = 1
    charts_y = 2
    pylab.figure(figsize=(charts_y, charts_x))
    pylab.subplots_adjust(hspace=.7)

    pylab.subplot(charts_y, charts_x, 1)
    pylab.title("Pixel Average")
    pylab.plot(averages)

    frequencies = scipy.fftpack.fftfreq(len(averages), d=1.0 / fps)

    pylab.subplot(charts_y, charts_x, 2)
    pylab.title("FFT")
    pylab.axis([0, 15, -2000000, 5000000])
    pylab.plot(frequencies, scipy.fftpack.fft(averages))

    pylab.show()

def display_head(set_x, set_y, n = 5):
    '''
    show some figures based on gray image matrixs
    
    @type set_x: TensorSharedVariable, 
    @param set_x: gray level value matrix of the
    
    @type set_y: TensorVariable, 
    @param set_y: label of the figures    
    
    @type n: int, 
    @param n: numbers of figure to be display, less than 10, default 5
    '''
    import pylab
    
    if n > 10: n = 10
    img_x = set_x.get_value()[0:n].reshape(n, 28, 28)
    img_y = set_y.eval()[0:n]
    
    for i in range(n): 
        pylab.subplot(1, n, i+1); 
        pylab.axis('off'); 
        pylab.title(' %d' % img_y[i])
        pylab.gray()
        pylab.imshow(img_x[i])

def plot_features(im, features, num_to_plot=100, colors=["blue"]):

  plt.imshow(im)

  for i in range(min(features.shape[0], num_to_plot)):
    x = features[i,0]
    y = features[i,1]
    scale = features[i,2]
    rot = features[i,3]


    color = colors[i % len(colors)]

    box = patches.Rectangle((-scale/2,-scale/2), scale, scale, 
      edgecolor=color, facecolor="none", lw=1)
    arrow = patches.Arrow(0, -scale/2, 0, scale, 
      width=10, edgecolor=color, facecolor="none")
    t_start = plt.gca().transData
    transform = mpl.transforms.Affine2D().rotate(rot).translate(x,y) + t_start

    box.set_transform(transform)
    arrow.set_transform(transform)
    plt.gca().add_artist(box)
    plt.gca().add_artist(arrow)


  plt.axis('off')
  plt.set_cmap('gray')
  plt.show()