def get_cov_ev(cov, plot=False):
    import numpy as np
    import pandas as pd
    cov = np.mat(cov)
    eig_vals, eig_vecs = np.linalg.eig(cov)
    df = pd.DataFrame()
    df['ev_x'] = eig_vecs[0, :].tolist()[0]
    df['ev_y'] = eig_vecs[1, :].tolist()[0]
    df['e_values'] = eig_vals
    df['e_values_sqrt'] = df.e_values.apply(np.sqrt)
    df.fillna(0, inplace=True)
    df = df.sort_values(by='e_values', ascending=False).reset_index()
    main_axis=df.ev_x[0] * df.e_values_sqrt[0], df.ev_y[0] * df.e_values_sqrt[0]
    angle = vector_to_angle(main_axis[0], main_axis[1])  # from -135 to 45 degrees
    angle = angle if angle >= -90 else 180+angle  # to make it from -90 to 90

    if plot:
        print('cov\n\t' + str(cov).replace('\n', '\n\t'))
        import pylab as plt
        samples = 100
        df_original_multi_samples = pd.DataFrame(np.random.multivariate_normal([0, 0], cov, samples), columns=['X', 'Y'])
        df_original_multi_samples.plot.scatter(x='X', y='Y', alpha=0.1)
        m = df_original_multi_samples.max().max()
        plt.axis([-m, m, -m, m])

        for i in [0, 1]:  # 0 is the main eigenvector
            if df.ev_x[i] * df.e_values_sqrt[i] + df.ev_y[i] * df.e_values_sqrt[i]:  # cannot plot 0 size vector
                plt.arrow(0, 0,
                          df.ev_x[i] * df.e_values_sqrt[i], df.ev_y[i] * df.e_values_sqrt[i],
                          head_width=0.15, head_length=0.15,
                          length_includes_head=True, fc='k', ec='k')
            else:
                print('zero length eigenvector, skipping')
        plt.show()
    return df, angle

 def arc(xx, yy, index, draw_dots=True):
     A = array([xx[0], yy[0]])
     B = array([xx[1], yy[1]])
     if draw_dots:
         plot(A[0], A[1], "ko")
         plot(B[0], B[1], "ko")
     AB = B-A
     AB_len = sqrt(sum(AB**2))
     AB = AB / AB_len
     pAB = array([AB[1], -AB[0]])
     AB_mid = (A+B)/2.
     if index == 1:
         R = AB_mid + pAB * AB_len/1.5
     else:
         R = AB_mid + pAB * AB_len/4.
     r = sqrt(sum((A-R)**2))
     P_arrow = R - pAB * r
     # Draw the arc from A to B centered at R
     phi1 = atan2((A-R)[1], (A-R)[0])
     phi2 = atan2((B-R)[1], (B-R)[0])
     n = 100 # number of divisions
     phi = linspace(phi1, phi2, n)
     xx = r*cos(phi)+R[0]
     yy = r*sin(phi)+R[1]
     plot(xx, yy, "k-", lw=2)
     x, x2 = xx[n/2-1:n/2+1]
     y, y2 = yy[n/2-1:n/2+1]
     dx = x2-x
     dy = y2-y
     arrow(x, y, dx/2, dy/2, fc="black", head_width=0.05)

def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300, optimize=True, plot=True):
    """
    Show an example of a multimodal error surface for Gaussian process
    regression. Gene 939 has bimodal behaviour where the noisy mode is
    higher.
    """

    # Contour over a range of length scales and signal/noise ratios.
    length_scales = np.linspace(0.1, 60., resolution)
    log_SNRs = np.linspace(-3., 4., resolution)

    try:import pods
    except ImportError:
        print 'pods unavailable, see https://github.com/sods/ods for example datasets'
        return
    data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
    # data['Y'] = data['Y'][0::2, :]
    # data['X'] = data['X'][0::2, :]

    data['Y'] = data['Y'] - np.mean(data['Y'])

    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF)
    if plot:
        pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
        ax = pb.gca()
        pb.xlabel('length scale')
        pb.ylabel('log_10 SNR')

        xlim = ax.get_xlim()
        ylim = ax.get_ylim()

    # Now run a few optimizations
    models = []
    optim_point_x = np.empty(2)
    optim_point_y = np.empty(2)
    np.random.seed(seed=seed)
    for i in range(0, model_restarts):
        # kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
        kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))

        m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern)
        m.likelihood.variance = np.random.uniform(1e-3, 1)
        optim_point_x[0] = m.rbf.lengthscale
        optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        # optimize
        if optimize:
            m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters)

        optim_point_x[1] = m.rbf.lengthscale
        optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        if plot:
            pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1] - optim_point_x[0], optim_point_y[1] - optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
        models.append(m)

    if plot:
        ax.set_xlim(xlim)
        ax.set_ylim(ylim)
    return m # (models, lls)

	def plot_profile(self) :
		self.logger.debug('Plotting Profile curves') 
		try :
			xnew = np.linspace(self.time[0],self.time[-1],num=50)
			#Function handle for generating plot interpolate points
			f1 = lambda(tck) : interpolate.splev(xnew,tck)
			ynew = np.array(map(f1,self.tcks)).T
			ax = pl.gca()
			color_list = ['r','g','b','c','y']
			ax.set_color_cycle(color_list)
			pl.plot(self.time,self.var,'x',xnew,ynew);
			# plot arrows 
			arrow_scale = 40.0
			dx = (self.time[-1]-self.time[0])/arrow_scale;
			dy = self.slopes.T * dx  # transpose operation here
			for ii in range(self.n_sample) : 
				t = self.time[ii]
				X = self.var[ii,:]
				DY = dy[:,ii]
				for jj in range(self.n_var) :
					pl.arrow(t,X[jj],dx,DY[jj],linewidth=1)
			pl.title('Plot of spline fitted biochem profile vs time')
			pl.xlabel('time')
			pl.ylabel('profile')
			pl.axis('tight')
			pl.show()
		except AttributeError :
			sys.stderr.writelines(
			"Define Profile before trying to plot..!")
			raise

	def add_point(event):
		x, y = event.xdata, event.ydata
		points.append((x, y))
		pylab.cla()
		pylab.scatter(*zip(*points))
		pylab.xlim(-10, 10)
		pylab.ylim(-10, 10)

		pylab.draw()
		if len(points) < 3: return

		c, r = three_point_circle(*points)
		cir = pylab.Circle(c, r)
		pylab.gca().add_patch(cir)

		for p in points:
			angle = angle_at_point(c, p)
			if angle < 0:
				angle += 2*np.pi
			if angle >= np.pi:
				angle = angle - np.pi
			print np.degrees(angle)
			dx, dy = np.array((np.cos(angle), np.sin(angle)))
			pylab.text(p[0], p[1], "%.2f"%np.degrees(angle))
			pylab.arrow(p[0], p[1], dx, dy)
		pylab.show()

def plotAssignment(idx, coords, centroids, warehouseCoord, output):
    V = max(idx) + 1
    cmap = plt.cm.get_cmap('Dark2')
    customerColors = [cmap(1.0 * idx[i] / V) for i in range(len(idx))]
    centroidColors = [cmap(1.0 * i / V) for i in range(V)]
    xy = coords

    pylab.scatter([warehouseCoord[0]], [warehouseCoord[1]])
    pylab.scatter(centroids[:,0], centroids[:,1], marker='o', s=500, linewidths=2, c='none')
    pylab.scatter(centroids[:,0], centroids[:,1], marker='x', s=500, linewidths=2, c=centroidColors)
    pylab.scatter(xy[:,0], xy[:,1], s=100, c=customerColors)

    for cluster in range(V):
        customerIndices = indices(cluster, idx)
        clusterCustomersCoords = [warehouseCoord] + [list(coords[i]) for i in customerIndices]

        N = len(clusterCustomersCoords)
        for i in range(N):
            j = (i+1) % N
            x = clusterCustomersCoords[i][0]
            y = clusterCustomersCoords[i][1]
            dx = clusterCustomersCoords[j][0] - clusterCustomersCoords[i][0]
            dy = clusterCustomersCoords[j][1] - clusterCustomersCoords[i][1]
            pylab.arrow(x, y, dx, dy, color=centroidColors[cluster], fc="k",
                        head_width=1.0, head_length=2.5, length_includes_head=True)

    pylab.savefig(output)
    pylab.close()

    def plot(self,plot_frame,**kwargs):
        import pylab as plb
        default_args = {'draw_frame':True,
                        'frame_head_width':20,
                        'contour_kwargs':{'edgecolor': 'none', 
                                          'linewidth': 0.0, 
                                          'facecolor': 'none'}}
        from pylab import plot,arrow
        lines = self.model.coords_from_frame(plot_frame)
        self.curves = list()
        #plot_args = kwargs.pop('plot_args',default_args)
        for line_key in lines.keys():
            try:
                element_args = kwargs['contour_kwargs'][line_key]
            except KeyError:
                
                element_args = default_args['contour_kwargs']
            line = lines[line_key]
            from matplotlib.patches import Polygon
            poly = Polygon(zip(line[0,:],line[1,:]),**element_args)
            #plb.plot([1,2,3,4])
            plb.gca().add_patch(poly,)

        if 'draw_frame' in kwargs.keys():
            if kwargs['draw_frame']:
                frame_args = dict()
                p = plot_frame['p']
                a1 = plot_frame['a1']
                a2 = plot_frame['a2']
                frame_args['color'] = 'g'
                frame_args['head_width'] = kwargs['frame_head_width']
                arrow(p[0],p[1],a1[0],a1[1],**frame_args)
                frame_args['color'] = 'b'
                frame_args['head_width'] = kwargs['frame_head_width']
                arrow(p[0],p[1],a2[0],a2[1],**frame_args)

def plot_ch():
    for job in jobs_orig:
        print "plane of", job.path
        pylab.clf()
        x_center = int((job.position(0)[0] + job.position(1)[0])/2)
        x_final = 50 + x_center
        #plane = np.concatenate((job.plane(y=50)[:, x_final:], 
        #                        job.plane(y=50)[:, :x_final]), axis=1)
        plane = job.plane(y=50)
        myplane = plane[plane < 0.0]
        p0 = myplane.min()
        p12 = np.median(myplane)
        p14 = np.median(myplane[myplane<p12])
        p34 = np.median(myplane[myplane>p12])
        p1 = myplane.max()
        contour_values = (p0, p14, p12, p34, p1)
        pylab.title(r'$u=%.4f,\  D=%.4f,\ Q=%i$ ' %
                ((job.u_x**2+job.u_y**2)**0.5, job.D_minus, job.ch_objects))
        car = pylab.imshow(plane, vmin=-0.001, vmax=0.0, 
                           interpolation='nearest')
        pylab.contour(plane, contour_values, linestyles='dashed', 
                                             colors='white')
        print job.u_x, job.u_y
        pylab.grid(True)
        pylab.colorbar(car)
        pylab.arrow(1, 1, 2*job.u_x/(job.u_x if job.u_x else 1.0), 2*job.u_y/(job.u_y if job.u_y else 1.0), width=0.5)
        #imgfilename = 'plane_r20-y50-u_x%.4fD%.4fch%03i.png' % \
        #              (job.u_x, job.D_minus, job.ch_objects)
        imgfilename = 'plane_%s.png' % job.job_id
        pylab.savefig(imgfilename)

def triangle_vector_space(S,T,r, steps = 21):
	"""
	This draws the phase space in term of the triangular description of phase space, given a game
	specified by S, T and r. Steps defines the number ofdecrete steps to take in each axis.
	
	Note
	====
	This was a mostly experimental idea, it is not currently in use

	"""
	##Initialise a population at a certain point, with small frequencies of mutants around it

	M = GP_evo3(S,T,r=r, steps = steps)

	pl.figure()
	pl.plot( [0,.5],[0,1], color = 'black', linewidth = 1 )
	pl.plot( [0.5,1],[1,0], color = 'black', linewidth = 1 )
	xa_s = np.linspace(0,1,101)
	pl.plot( xa_s, map(M.phi_R,xa_s), '--', color = 'black', linewidth = 2.5 )
	pl.xlabel( '$x_A$', fontsize = 30 )
	pl.ylabel( '$\\varphi$', fontsize = 30 )
	for i in xrange(steps):
		for j in xrange(steps):

			arrow = delta_x(i,j,M)
			#print arrow
			pl.arrow( *arrow )

	pl.show()

def draw_bar(xs, y, dat):
        global scale

        for ix in range(len(xs)):
                if dat[ix] > 0: color = "red"
                else: color = "blue"
                pl.arrow(xs[ix], y/2., scale*dat[ix], 0, color=color)

def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, optim_iters=300):
    """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""

    # Contour over a range of length scales and signal/noise ratios.
    length_scales = np.linspace(0.1, 60.0, resolution)
    log_SNRs = np.linspace(-3.0, 4.0, resolution)

    data = GPy.util.datasets.della_gatta_TRP63_gene_expression(gene_number)
    # data['Y'] = data['Y'][0::2, :]
    # data['X'] = data['X'][0::2, :]

    data["Y"] = data["Y"] - np.mean(data["Y"])

    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.rbf)
    pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
    ax = pb.gca()
    pb.xlabel("length scale")
    pb.ylabel("log_10 SNR")

    xlim = ax.get_xlim()
    ylim = ax.get_ylim()

    # Now run a few optimizations
    models = []
    optim_point_x = np.empty(2)
    optim_point_y = np.empty(2)
    np.random.seed(seed=seed)
    for i in range(0, model_restarts):
        # kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
        kern = GPy.kern.rbf(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))

        m = GPy.models.GPRegression(data["X"], data["Y"], kernel=kern)
        m["noise_variance"] = np.random.uniform(1e-3, 1)
        optim_point_x[0] = m["rbf_lengthscale"]
        optim_point_y[0] = np.log10(m["rbf_variance"]) - np.log10(m["noise_variance"])

        # optimize
        m.optimize("scg", xtol=1e-6, ftol=1e-6, max_f_eval=optim_iters)

        optim_point_x[1] = m["rbf_lengthscale"]
        optim_point_y[1] = np.log10(m["rbf_variance"]) - np.log10(m["noise_variance"])

        pb.arrow(
            optim_point_x[0],
            optim_point_y[0],
            optim_point_x[1] - optim_point_x[0],
            optim_point_y[1] - optim_point_y[0],
            label=str(i),
            head_length=1,
            head_width=0.5,
            fc="k",
            ec="k",
        )
        models.append(m)

    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    return m  # (models, lls)

	def scatter_particles(self, state):
		for particle in state["particles"]:
			x, y = particle["x"]
			vx, vy = particle["v"]
			if particle["id"] == state["best_id"]:
				pl.scatter(x, y, c='y', s=35)
			else:
				pl.scatter(particle["x"][0], particle["x"][1], c='k')
			pl.arrow(x, y, vx, vy, shape="full", head_width=0.03, width=0.00001, alpha=0.3)

def test_optimal_angle_filter():
	image = pylab.imread(sys.argv[1])
	image = np.mean(image, axis=2)
	image = image[::-1]
	pylab.gray()

	f = AngleFilter((7, 7))

	components = f.get_component_values(image, mode='same')

	angles = np.radians(range(0, 180, 1))
	def best_angle_value(c):
		values = f.basis.response_at_angle(c, angles)
		return angles[np.argmax(values)], np.max(values)

	for y, x in np.ndindex(components.shape[:2]):
		if y%4 != 0: continue
		if x%4 != 0: continue
		maxval = f.basis.max_response_value(components[y,x])
		if maxval < 2: continue
		maxang_an = f.basis.max_response_angle(components[y,x])
		maxang, maxval = best_angle_value(components[y,x])
		pylab.scatter(x, y)
		dy = -5.0
		dx, dy = np.array((np.cos(maxang), np.sin(maxang)))*10
		#dx = np.tan(maxang_an)*dy
		#d = d/np.linalg.norm(d)*3
		pylab.arrow(x, y, dx, dy, color='blue')

		#d = np.array((-np.sin(maxang_an), -np.cos(maxang_an)))
		#d = d/np.linalg.norm(d)*3
		#pylab.arrow(x, y, d[0], d[1], color='green')

	
		#pylab.plot(x, y, '.')

	pylab.imshow(image, origin="lower")
	#pylab.xlim(0, components.shape[1])
	#pylab.ylim(components.shape[0], 0)
	pylab.show()
	return
	
	#pylab.subplot(1,2,1)
	#pylab.imshow(image)
	#pylab.subplot(1,2,2)
	filtered = np.zeros(image.shape)
	for y, x in np.ndindex(components.shape[:2]):
		maxval = f.basis.max_response_value(components[y,x])
		minval = f.basis.min_response_value(components[y,x])
		#print maxval, minval
		filtered[y,x] = maxval
		#filtered[y, x] = best_angle_value(components[y,x])
	#pylab.hist(filtered.flatten())
	pylab.imshow(filtered > 3)
	pylab.show()

def imshow_cube_image(image, header=None, compass=True, blackhole=True, 
                      cmap=None):
    """
    Call imshow() to plot a cube image. Make sure it is already 
    masked. Also pass in the header to do the compass rose calculations.
    """
    # Setup axis info (x is the 2nd axis)
    xcube = (np.arange(image.shape[1]) - m31pos[0]) * 0.05
    ycube = (np.arange(image.shape[0]) - m31pos[1]) * 0.05
    xtickLoc = py.MultipleLocator(0.4)

    if cmap is None:
        cmap = py.cm.jet

    # Plot the image.
    py.imshow(image, 
              extent=[xcube[0], xcube[-1], ycube[0], ycube[-1]],
              cmap=cmap)
    py.gca().get_xaxis().set_major_locator(xtickLoc)
    
    py.xlabel('X (arcsec)')
    py.ylabel('Y (arcsec)')

    # Get the spectrograph position angle for compass rose
    if compass is True:
        if header is None:
            pa = paSpec
        else:
            pa = header['PA_SPEC']


        # Make a compass rose
        cosSin = np.array([ math.cos(math.radians(pa)), 
                            math.sin(math.radians(pa)) ])
        arr_base = np.array([ xcube[-1]-0.5, ycube[-1]-0.6 ])
        arr_n = cosSin * 0.2
        arr_w = cosSin[::-1] * 0.2
        py.arrow(arr_base[0], arr_base[1], arr_n[0], arr_n[1],
                 edgecolor='w', facecolor='w', width=0.03, head_width=0.08)
        py.arrow(arr_base[0], arr_base[1], -arr_w[0], arr_w[1],
                 edgecolor='w', facecolor='w', width=0.03, head_width=0.08)
        py.text(arr_base[0]+arr_n[0]+0.1, arr_base[1]+arr_n[1]+0.1, 'N', 
                color='white', 
                horizontalalignment='left', verticalalignment='bottom')
        py.text(arr_base[0]-arr_w[0]-0.15, arr_base[1]+arr_w[1]+0.1, 'E', 
                color='white',
                horizontalalignment='right', verticalalignment='center')

    if blackhole is True:
        py.plot([0], [0], 'ko')


    py.xlim([-0.5, 0.6])
    py.ylim([-2.0, 1.8])

def drawScaledEigenvectors(X,Y, eigVectors, eigVals, theColor='k'):
    """ Draw scaled eigenvectors starting at (X,Y)"""
    
    # For each eigenvector
    for col in xrange(eigVectors.shape[1]):
        
        # Draw it from (X,Y) to the eigenvector length
        # scaled by the respective eigenvalue
        pylab.arrow(X,Y,eigVectors[0,col]*eigVals[col],
                    eigVectors[1,col]*eigVals[col],
                    width=0.01, color=theColor)

def arrow(start=(0, 0), vel=(0, (0, 0)), color='b'):
    if vel == (0, (0, 0)):
        xr, yr = 0.1, 0.1
    else:
        if vel[0] == 0:
            sp = 0.1
        else:
            sp = vel[0]
        x = vel[1][0]
        y = vel[1][1]
        xr = sp * x * qq(x, y)
        yr = sp * y * qq(x, y)
    pylab.arrow(start[0], start[1], xr, yr, head_width=0.05, head_length=0.1, color=color)

def show_angle(angle,power):
    ARROW_STEP = 40
    kernel = np.ones((ARROW_STEP,ARROW_STEP))
    rowFilter = gaussian(ARROW_STEP,ARROW_STEP/5)
    colFilter = rowFilter
    
    gb180 = cmt.get_cmap('gb180')
    
    #X = np.arange(0,angle.shape(-1),ARROW_STEP)
    #Y = np.arange(0,angle.shape(-2),ARROW_STEP)

    x = np.matrix(np.arange(ARROW_STEP/2,angle.shape[-1],ARROW_STEP))
    y = np.transpose(np.matrix(np.arange(ARROW_STEP/2,angle.shape[-2],ARROW_STEP)))

    X = np.array((0*y+1)*x)
    Y = np.array(y*(0*x+1))
    
    #u = convolve2d(sin(angle),kernel,mode='same')
    #v = convolve2d(cos(angle),kernel,mode='same')
    #p = convolve2d(power,kernel,mode='same')
    
    u = sepfir2d(np.sin(angle),rowFilter,colFilter)
    v = sepfir2d(np.cos(angle),rowFilter,colFilter)
    p = sepfir2d(power,rowFilter,colFilter)
    
    U = u[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP]*p[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP]
    V = v[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP]*p[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP]
    
    #U = sin(angle[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP])*(power[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP])
    #V = cos(angle[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP])*(power[ARROW_STEP/2::ARROW_STEP,ARROW_STEP/2::ARROW_STEP])

    #X = X[(power[::ARROW_STEP,::ARROW_STEP]>.016)]
    #Y = Y[(power[::ARROW_STEP,::ARROW_STEP]>.016)]
    #U = U[(power[::ARROW_STEP,::ARROW_STEP]>.016)]
    #V = V[(power[::ARROW_STEP,::ARROW_STEP]>.016)]
    
    ua = ARROW_STEP/1.5*np.nansum(np.sin(angle)*power)/np.nansum(power)
    va = -ARROW_STEP/1.5*np.nansum(np.cos(angle)*power)/np.nansum(power)
    xc, yc = angle.shape[-1]/2, angle.shape[-2]/2
    
    pylab.imshow(angle,cmap=gb180)
    #pylab.imshow(angle,cmap='hsv')
    pylab.quiver(X,Y,U,V,pivot='middle',color='w',headwidth=1,headlength=0)
    pylab.arrow(xc,yc,ua,va,color='w',linewidth=2)
    pylab.arrow(xc,yc,-ua,-va,color='w',linewidth=2)
    ax=pylab.gca()
    ax.set_axis_off()
    pylab.show()

    A = np.arctan2(-ua,va)
    return A

def riemann():

    
    # grid info
    xmin = 0.0
    xmax = 1.0

    nzones = 1
    ng = 0
    
    gr = gpu.grid(nzones, xmin=xmin, xmax=xmax)


    #------------------------------------------------------------------------
    # plot a domain without ghostcells
    gpu.drawGrid(gr)

    gpu.labelCenter(gr, 0, r"$i$")

    gpu.labelCellCenter(gr, 0, r"$q_i$")
    
    gpu.markCellLeftState(gr, 0, r"$q_{i-1/2,R}^{n+1/2}$", color="r")
    gpu.markCellRightState(gr, 0, r"$q_{i+1/2,L}^{n+1/2}$", color="r")
    

    pylab.arrow(gr.xc[0]-0.05*gr.dx, 0.5, -0.13*gr.dx, 0, 
                shape='full', head_width=0.075, head_length=0.05, 
                lw=1, width=0.01,
                edgecolor="none", facecolor="r",
                length_includes_head=True, zorder=100)

    pylab.arrow(gr.xc[0]+0.05*gr.dx, 0.5, 0.13*gr.dx, 0, 
                shape='full', head_width=0.075, head_length=0.05, 
                lw=1, width=0.01,
                edgecolor="none", facecolor="r",
                length_includes_head=True, zorder=100)
    

    pylab.xlim(gr.xl[0]-0.25*gr.dx,gr.xr[2*ng+nzones-1]+0.25*gr.dx)
    # pylab.ylim(-0.25, 0.75)
    pylab.axis("off")

    pylab.subplots_adjust(left=0.05,right=0.95,bottom=0.05,top=0.95)

    f = pylab.gcf()
    f.set_size_inches(4.0,2.5)


    pylab.savefig("states.png")
    pylab.savefig("states.eps")

def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000):
    """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""

    # Contour over a range of length scales and signal/noise ratios.
    length_scales = np.linspace(0.1, 60., resolution)
    log_SNRs = np.linspace(-3., 4., resolution)

    data = GPy.util.datasets.della_gatta_TRP63_gene_expression(gene_number)
    # Sub sample the data to ensure multiple optima
    #data['Y'] = data['Y'][0::2, :]
    #data['X'] = data['X'][0::2, :]

    # Remove the mean (no bias kernel to ensure signal/noise is in RBF/white)
    data['Y'] = data['Y'] - np.mean(data['Y'])

    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.rbf)
    pb.contour(length_scales, log_SNRs, np.exp(lls), 20)
    ax = pb.gca()
    pb.xlabel('length scale')
    pb.ylabel('log_10 SNR')

    xlim = ax.get_xlim()
    ylim = ax.get_ylim()

    # Now run a few optimizations
    models = []
    optim_point_x = np.empty(2)
    optim_point_y = np.empty(2)
    np.random.seed(seed=seed)
    for i in range(0, model_restarts):
        kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + GPy.kern.white(1,variance=np.random.exponential(1.))

        m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern)
        optim_point_x[0] = m.get('rbf_lengthscale')
        optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));

        # optimize
        m.ensure_default_constraints()
        m.optimize(xtol=1e-6,ftol=1e-6)

        optim_point_x[1] = m.get('rbf_lengthscale')
        optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));

        pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1]-optim_point_x[0], optim_point_y[1]-optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
        models.append(m)

    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    return (models, lls)

def snapshot(t, pos, vel, colors, arrow_scale=0.2):
    global img
    pylab.cla()
    pylab.axis([0, 1, 0, 1])
    pylab.setp(pylab.gca(), xticks=[0, 1], yticks=[0, 1])
    for (x, y), (dx, dy), c in zip(pos, vel, colors):  # putting position, velocity, and colors arrays into 3D tuple.
        dx *= arrow_scale  # rescale vel_x by multiplying w/ arrow_scale
        dy *= arrow_scale  # rescale vel_y by multiplying w/ arrow_scale
        circle = pylab.Circle((x, y), radius=sigma, fc=c)  # fc sets color of circle
        pylab.gca().add_patch(circle)
    pylab.arrow(x, y, dx, dy, fc="k", ec="k", head_width=0.05, head_length=0.05)
    pylab.text(
        0.5, 1.03, "t = %.2f" % t, ha="center"
    )  # .5 for middle of box on x-axis, 1.03 is just above box on y-axis
    pylab.savefig(os.path.join(output_dir, "%d.png" % img))
    img += 1