def work_on_area(t, do_plot=True):
    if do_plot:
        plt.figure(figsize=(4,4))
    lst = []
    for i in range(N):
        for j in range(N):
            ff[i,j] = np.std(ens([s[t], h[t]], [k_arg[i],k_arg[j]]) - m_fc[t])
            lst.append([k_arg[i],k_arg[j],ff[i,j]])        
    lst = np.array(lst)
    ff_min = np.min(ff)
    ff_sel = lst[lst[:,2].flatten() < ff_min*1.1][:,0:2].transpose()
    min_x = k_arg[np.argmin(ff) // ff.shape[1]]
    min_y = k_arg[np.argmin(ff) % ff.shape[1]]
    pp.fit(ff_sel.transpose())
    if do_plot:
        plt.contour(k_arg, k_arg, ff.transpose() - np.min(ff))
        plt.plot(ff_sel[0], ff_sel[1], 'o')
        plt.plot(min_x, min_y, 'or')    
        plt.quiver([min_x, min_x], [min_y, min_y],
                   pp.components_[:, 0].flatten() * pp.explained_variance_,
                   pp.components_[:, 1].flatten() * pp.explained_variance_,
                   scale=1.0)
        plt.xlim(lims)
        plt.ylim(lims)
    return np.concatenate(([min_x, min_y],pp.components_.flatten(),pp.explained_variance_))

	def __call__(self,u,w,bx,by,bz,b2,t):
		q = 8

		map = cm.red_blue()
		if self.x == None:
			nx = u.shape[2]
			nz = u.shape[0]
			self.x,self.y = np.meshgrid(range(nx),range(nz))
		x,y = self.x,self.y

		avgu = np.average(u,1)
		avgw = np.average(w,1)
		avgbx = np.average(bx,1)
		avgby = np.average(by,1)
		avgbz = np.average(bz,1)
		avgb2 = np.average(b2,1)
		avgt = np.average(t,1)

		plt.subplot(121)
		plt.imshow(avgt,cmap=map,origin='lower')
		plt.colorbar()
		plt.quiver(x[::q,::q],y[::q,::q],avgu[::q,::q],avgw[::q,::q])
		plt.title('Tracer-Vel')
		plt.axis("tight")

		plt.subplot(122)
		plt.imshow(avgby,cmap=map,origin='lower')
		plt.colorbar()
		plt.quiver(x[::q,::q],y[::q,::q],avgbx[::q,::q],avgbz[::q,::q])
		plt.title('By-Twist')
		plt.axis("tight")

def _decorate_contour_segment(data, stride=1, options={}, tomax=True, labelled=False, outline=None, aspect=1):
	default_options = {'scale': 0.2,
			'scale_units': 'dots',
			'headaxislength': 2,
			'headlength': 2,
			'headwidth': 2,
			'minshaft': 1,
			'units': 'dots',
			#'angles': 'xy',
			'edgecolor': outline,
			'linewidth': 0 if outline is None else 0.2
		}
	default_options.update(options)

	x = data[::stride,0]
	y = data[::stride,1]

	sign = 1 if tomax else -1
	dx = -sign*np.diff(y)*aspect
	dy = sign*np.diff(x)
	l = np.sqrt(dx**2+dy**2)
	dx /= l
	dy /= l

	x = 0.5*(x+np.roll(x,-1))
	y = 0.5*(y+np.roll(y,-1))

	if labelled:
		x,y,dx,dy = x[1:-2], y[1:-2], dx[1:-1], dy[1:-1]
	else:
		x,y = x[:-1], y[:-1]

	plt.quiver(x, y, dx, dy, **default_options)

def quiverPlot(psi, orient, saveDest):
    """ Generates a quiver plot that shows channel orientation
    and singularity index response strength.

    Inputs:
    psi -- singularity index response
    orient -- local orientation at each spatial location (x,y)
    saveDest -- output figure save destination

    Returns:
    None (saves the figure at saveDest)
    """

    # downsample
    psi_s = psi[::4,::4]
    orient_s = orient[::4,::4]

    U = -psi_s * np.sin(orient_s)
    V = psi_s * np.cos(orient_s)

    R, C = psi.shape
    aspect_ratio = float(R)/C

    plt.figure(figsize=(10, 10*aspect_ratio))
    ax = plt.gca()
    ax.invert_yaxis()
    plt.quiver(U, V, scale=4, width=0.001, pivot='mid', headwidth=0, minlength=0)
    plt.axis('off')
    plt.savefig(saveDest)

def test_complete():
    fig = plt.figure('Figure with a label?', figsize=(10, 6))

    plt.suptitle('Can you fit any more in a figure?')

    # make some arbitrary data
    x, y = np.arange(8), np.arange(10)
    data = u = v = np.linspace(0, 10, 80).reshape(10, 8)
    v = np.sin(v * -0.6)

    plt.subplot(3, 3, 1)
    plt.plot(list(xrange(10)))

    plt.subplot(3, 3, 2)
    plt.contourf(data, hatches=['//', 'ooo'])
    plt.colorbar()

    plt.subplot(3, 3, 3)
    plt.pcolormesh(data)

    plt.subplot(3, 3, 4)
    plt.imshow(data)

    plt.subplot(3, 3, 5)
    plt.pcolor(data)

    plt.subplot(3, 3, 6)
    plt.streamplot(x, y, u, v)

    plt.subplot(3, 3, 7)
    plt.quiver(x, y, u, v)

    plt.subplot(3, 3, 8)
    plt.scatter(x, x**2, label='$x^2$')
    plt.legend(loc='upper left')

    plt.subplot(3, 3, 9)
    plt.errorbar(x, x * -0.5, xerr=0.2, yerr=0.4)

    ###### plotting is done, now test its pickle-ability #########

    # Uncomment to debug any unpicklable objects. This is slow (~200 seconds).
#    recursive_pickle(fig)

    result_fh = BytesIO()
    pickle.dump(fig, result_fh, pickle.HIGHEST_PROTOCOL)

    plt.close('all')

    # make doubly sure that there are no figures left
    assert_equal(plt._pylab_helpers.Gcf.figs, {})

    # wind back the fh and load in the figure
    result_fh.seek(0)
    fig = pickle.load(result_fh)

    # make sure there is now a figure manager
    assert_not_equal(plt._pylab_helpers.Gcf.figs, {})

    assert_equal(fig.get_label(), 'Figure with a label?')

 def test_dir_segment_complet(self):
     plt.figure('test segment complet')
     X, Y = np.mgrid[-20:20:20j, -20:20:20j]
     U1, V1 = vfl.dir_segment(X, Y, -5, 0, 5, 0)
     U2, V2 = vfl.dir_segment_extremity(X, Y, -5, 0, 5, 0)
     plt.quiver(X, Y, U1 + U2, V1 + V2)
     plt.show()

def three_d_contour(file,show):

    xloc = raw_input("xlocation(two_d_initial_x_value.txt):") or "two_d_initial_x_value.txt"
    yloc = raw_input("ylocation(two_d_initial_y_value.txt):") or "two_d_initial_y_value.txt"
    uloc = raw_input("ulocation(two_d_cavity_u_results.txt):") or "two_d_cavity_u_results.txt"
    vloc = raw_input("vlocation(two_d_cavity_v_results.txt):") or "two_d_cavity_v_results.txt"
    ploc = raw_input("plocation(two_d_cavity_p_results.txt):") or "two_d_cavity_p_results.txt"
    u=np.loadtxt(uloc)
    v=np.loadtxt(vloc)
    p=np.loadtxt(ploc)
    X=np.loadtxt(xloc)
    Y=np.loadtxt(yloc)
    X,Y=np.meshgrid(X,Y)
#    print X,Y
    fig = plt.figure()
    plt.contourf(X,Y,p,alpha=0.5)    ###plnttong the pressure field as a contour
    plt.colorbar()
    plt.contour(X,Y,p)               ###plotting the pressure field outlines
    plt.quiver(X[::3,::3],Y[::3,::3],u[::3,::3],v[::3,::3]) ##plotting velocity
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title(file)
    sep = '.'
    filename = file.split(sep, 1)[0]
    plt.savefig(filename+'.png')   # save figure as .png file
    if show:
        # print('test')
        plt.show()
    plt.close(fig)

def contour_plot(ax, isub, var, U, V, x_range, y_range, levels_var=None, annot=None):
    X, Y = np.meshgrid(x_range, y_range)
    plt.subplots_adjust(hspace=0.0)
    plt.subplots_adjust(wspace=0)
    temp = ticker.MaxNLocator(3)
    # if isub < 2:
    #    ax.set_xticklabels(())
    if isub % 2 != 0:
        ax.set_yticklabels(())
    # contour plot for v field
    CS = plt.contourf(X, Y, var, cmap=plt.cm.Blues, alpha=0.7, levels=levels_var)
    # plotting vectors from U and V fields
    if (U ** 2 + V ** 2).max() > 0:
        plt.quiver(X[::20, ::20], Y[::20, ::20], U[::20, ::20], V[::20, ::20], scale=15)
    # cbar = plt.colorbar(CS, fraction=0.04, format="%.1e")
    CS_l = plt.contour(CS, levels=CS.levels, colors="k", linewidths=0.2)

    ax.set_ylim(-8.5, 8.5)
    ax.set_xlim(-8.5, 8.5)
    ax.annotate(
        annot,
        xy=(0.01, 0.97),
        xycoords="axes fraction",
        fontsize=12,
        horizontalalignment="left",
        verticalalignment="top",
    )

def display_velocity( q , p ,mu , mu2 = None):
 W = 5*SIGMA
 res = 30
 N_nodes = res**2
 store = np.outer( np.linspace(-W,W , res), np.ones(res) )
 nodes = np.zeros( [N_nodes , DIM] )
 nodes[:,0] = np.reshape( store , N_nodes )
 nodes[:,1] = np.reshape( store.T , N_nodes )
 K,DK,D2K,D3K = jpf.derivatives_of_kernel(nodes , q)
 vel_field = np.einsum('ijab,jb->ia',K,p) - np.einsum('ijabc,jbc->ia',DK,mu)
 if mu2 != None:
  vel_field = vel_field + np.einsum('ijabcd,jbcd->ia',D2K,mu2)

 U = vel_field[:,0]
 V = vel_field[:,1]
 plt.figure()
 plt.quiver( nodes[:,0] , nodes[:,1] , U , V , scale=10 )
 plt.plot(q[:,0],q[:,1],'ro')
 for i in range(0,N):
     if np.linalg.norm(p[i]) < 1e-4: continue
     plt.arrow(q[i,0], q[i,1], p[i,0], p[i,1], head_width=0.1, head_length=0.2, lw = 4.0, fc='b', ec='b')
     plt.arrow(q[i,0], q[i,1], p[i,0], p[i,1], head_width=0.1, head_length=0.2, lw = 2.0, fc='w', ec='w')
 plt.axis('equal')
 plt.axis([- W, W,- W, W ])
 return plt.gcf()

def plotArrows( 
	# not yet tested
	prefix, dX, dY, X, Y,  slices, BG=None, C=None, 
	extent=None, zs = None, by=2,  figsize=default_figsize, 
	cmap=default_cmap, interpolation=default_interpolation, vmin=None, vmax=None, cbar=False,
	atoms=None, bonds=None, atomSize=default_atom_size
	):
	for ii,i in enumerate(slices):
		print " plotting ", i
		plt.figure( figsize=figsize )
		#plt.plt.quiver( dX, dY, X, Y, C=C, width=width, scale=scale )
		plt.quiver(  Xs[::by,::by], Ys[::by,::by],  dX[::by,::by], dY[::by,::by], color = 'k', headlength=10, headwidth=10, scale=15 )
		if BG is not None:
			plt.imshow    ( BG[i,:,:], origin='image',  interpolation=interpolation, cmap=cmap, extent=extent, vmin=vmin, vmax=vmax  )
			if cbar:
				plt.colorbar()
		plotGeom( atoms, bonds, atomSize=atomSize )
		plt.xlabel(r' Tip_x $\AA$')
		plt.ylabel(r' Tip_y $\AA$')
		if zs is None:
			plt.title( r"iz = %i" %i  )
		else:
			plt.title( r"Tip_z = %2.2f $\AA$" %zs[i]  )
		plt.savefig( prefix+'_%3.3i.png' %i, bbox_inches='tight' )
		plt.close()

def plot_vector_grid(filename, x_grid, y_grid):
    try:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
    except ImportError:
        print 'Plotting diasbled'
        return

    fig = plt.figure()
    x_coords = []
    y_coords = []
    z_coords = []
    for itt_x in range(len(x_grid[0])):
        for itt_y in range(len(x_grid)):
            x_coords.append(itt_x)
            y_coords.append(itt_y)
            z_coords.append((
                x_grid[itt_y][itt_x],
                y_grid[itt_y][itt_x]
            ))

    plt.quiver(x_coords,
               y_coords,
               [x[0] for x in z_coords],
               [y[1] for y in z_coords])
    fig.savefig(filename)

def weight_quiver(weights, color='c'):
    plt.quiver(weights[0, :-1],
               weights[1, :-1],
               weights[0, 1:] - weights[0, :-1],
               weights[1, 1:] - weights[1, :-1],
               scale_units='xy', angles='xy', scale=1,
               color=color)

def testLinVortList():
	"""Test function for a linear vortex sheet list"""
	V=linVortList()
	Npoints=50
	lowLim=numpy.array([0.,0.])
	upLim=numpy.array([1.,1.])
	for i in range(Npoints):
		p1=(lowLim+(float(i)/Npoints*(upLim-lowLim)))
		p2=(lowLim+(float(i+1)/Npoints*(upLim-lowLim)))
		V.addLinVortex(0.1,0.1,p1,p2)
	
	X,Y = numpy.meshgrid(numpy.arange(-2.,2.,0.2),numpy.arange(-2,2,0.2))
	u=numpy.copy(X)
	v=numpy.copy(X)
	
	for i in range(len(X.flatten())):
		vel=V.fieldEffect([X.flatten()[i],Y.flatten()[i]])
		u.ravel()[i]=vel[0]
		v.ravel()[i]=vel[1]
	
	plt.figure()
	plt.quiver(X,Y,u,v)
	plt.title('Vector field due to multiple linear vortex sheets')
	plt.plot([lowLim[0],upLim[0]],[lowLim[1],upLim[1]])
	plt.show()
	return()

def testVortexList(strength=1.0,blobType=0,delta=0.0):
	"""Test function for a point vortex list"""
	V=vortexList()
	Npoints=50
	lowLim=numpy.array([0.,0.])
	upLim=numpy.array([1.,1.])
	for i in range(Npoints):
		[x,y]=(lowLim+(float(i)/Npoints*(upLim-lowLim)))
		V.addVortex([x,y], strength,blobType,delta)
	
	X,Y = numpy.meshgrid(numpy.arange(-2.,2.,0.2),numpy.arange(-2,2,0.2))
	u=numpy.copy(X)
	v=numpy.copy(X)
	
	for i in range(len(X.flatten())):
		vel=V.fieldEffect([X.flatten()[i],Y.flatten()[i]])
		u.ravel()[i]=vel[0]
		v.ravel()[i]=vel[1]
	
	plt.figure()
	plt.quiver(X,Y,u,v)
	plt.title('Vector field due to multiple vortices')
	plt.scatter([lowLim[0],upLim[0]],[lowLim[1],upLim[1]])
	plt.show()
	return()

def plot_arrow(vectordata,imname,imdim):
    """ plot arrow plot """

    # plot slope direction
    rows = imdim[0]
    cols = imdim[1]
    fig = plt.figure(figsize=(rows/12.0,cols/12.0))
    fig.patch.set_alpha(0.0)
    fig.subplots_adjust(left=0.0,bottom=0.0,top=1.0,right=1.0)
    [X,Y,U,V] = vectordata
    if rows > 50 or cols > 50:
        #print "generate simple arrow plot"
        skip_p = int((rows*cols)/1000.0) 
        plt.quiver(X[::skip_p], Y[::skip_p], U[::skip_p], V[::skip_p],color='w')
    else:
        plt.quiver(X, Y, U, V,color='w')

    #plt.plot(X,Y,"bo")
    plt.axis('off')
    xup = np.unique(X)
    xstep = abs(xup[0] - xup[1])
    yup = np.unique(Y)
    ystep = abs(yup[0] - yup[1])
    plot_range = [np.min(X)-xstep,np.max(X)+xstep,np.min(Y)-ystep,np.max(Y)+ystep]
    plt.axis(plot_range)
    plt.savefig(imname + ".png", format="PNG",aspect="auto",transparent=False) # keep the white background

    plt.close(fig)

    return plot_range

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()
    plt.quiver(Q[0]-Q[1], Q[3]-Q[2])
    plt.show()

 def test_ligne(self):
     plt.figure('LIGNE ATTRACTIVE')
     X, Y = np.mgrid[-100:100:40j, -100:100:40j]
     U, V = vfl.ligne(
         X, Y, 0, 0, xb=50, yb=0, K=0.2, R=20, effect_range=40)
     plt.quiver(X, Y, U, V, scale=2)
     plt.show()

def plot_xwt_wavetransf(power, time, wa, T, S, sig95, pangle, time_base, scalemin=0, scalemax=6, ylabel='Pressure (mb)', plot_percentile=False):
    """plotting WaveTransform Power with confidence interval contour and phase vectors"""

    fig = plt.figure(10)
    ax = plt.subplot(1,1,1)
    if plot_percentile:
        #use following to contour at "percentiles variances" when using non-normalized data to match web output
        csf =plt.contourf(T, S, power, levels=[ 0, stats.scoreatpercentile(power, 25), stats.scoreatpercentile(power, 50),
                                           stats.scoreatpercentile(power, 75), stats.scoreatpercentile(power, 95), 
                                           stats.scoreatpercentile(power, 100)], colors=bmap)
    else:
        #use following to contour at "normalized variances" BAMS
        csf =plt.contourf(T, S, power, levels=[ 0, .2,.4,.6,.8,1], colors=bmap)
    cbar = plt.colorbar(pad=.1, shrink=.5, format='%.4f', extend='both') #move and shrink colorbar
    levels = [-99, 1] # values greater than 1 are significant
    plt.contour(T, S, sig95,levels, colors='black', linewidths=1)
    ax.set_yscale('log')
    ax.grid(True)
    
    # plot phase relationship 
    arr_dens = [60, 30]
    arr_densx = np.round( len(time) / arr_dens[0] )
    arr_densy = np.round( len(wa.scales) / arr_dens[1] )
    if arr_dens == 0:
        arr_dens = 1
        
    plt.quiver(T[::arr_densy,::arr_densx],S[::arr_densy,::arr_densx],(np.cos(pangle))[::arr_densy,::arr_densx],(np.sin(pangle))[::arr_densy,::arr_densx],\
    width=.00125, headwidth=4, headlength=4, alpha=0.6, color='k')

    # put the ticks at powers of 2 in the scale
    ticks = np.unique(2 ** np.floor(np.log2(wa.scales)))[1:]
    ax.yaxis.set_ticks(ticks)
    ax.yaxis.set_ticklabels(ticks.astype(str))
    ax.set_ylim(scalemax, scalemin)
    ax.set_ylabel('scales')

    # second y scale with equivalent fourier periods to scales
    # except with the ticks at the powers of 2
    ax_fourier = ax.twinx()
    ax_fourier.set_yscale('log')
    # match the fourier ticks to the scale ticks
    ax_fourier.set_yticks(ticks)
    ax_fourier.set_yticklabels(ticks.astype(str))
    ax_fourier.set_ylabel('fourier period (%s)' % time_base )
    fourier_lim = [wa.fourier_period(i) for i in ax.get_ylim()]
    ax_fourier.set_ylim(fourier_lim)

    ax.xaxis.set_major_formatter(DateFormatter('%Y-%m-%d'))
    fig.autofmt_xdate()

    # shade the region between the edge and coi
    C, S = wa.coi
    ax.fill_between(x=C, y1=S, y2=wa.scales.max(), color='gray', alpha=0.5)
    ax.set_xlim(time.min(), time.max())
        
    #plt.show()
    DefaultSize = fig.get_size_inches()
    fig.set_size_inches( (DefaultSize[0]*2, DefaultSize[1]) )
    
    return (plt, fig)

def plot_field(X, Y, U, V, filename):
    '''
    Function to plot the potential field.

    Args:
        X (numpy.ndarray): X component of the sample points.
        Y (numpy.ndarray): Y component of the sample points.
        U (numpy.ndarray): X component of field at sample points.
        V (numpy.ndarray): Y component of field at sample points.
    '''
    # Generate plot.
    padding = 0.5
    plt.figure()
    plt.quiver(X, Y, U, V,
               color='#007ce8',
               units='x',
               pivot='tail')
    plt.axis('equal')
    plt.axis([np.amin(X) - padding,
              np.amax(X) + padding,
              np.amin(Y) - padding,
              np.amax(Y) + padding])
    # plt.savefig("potential_field_back1.svg", format='svg')
    plt.savefig(filename + ".svg", format='svg')
    plt.show()

def render_still_vector_grid2(filename_x, filename_y, output_filename, **kwargs):
    has_xtick = True
    has_ytick = True
    if 'has_xtick' in kwargs:
        has_xtick = kwargs['has_xtick']
    if 'has_ytick' in kwargs:
        has_ytick = kwargs['has_ytick']
    grid_data_u = np.load(filename_x)
    grid_data_v = np.load(filename_y)
    nx = len(grid_data_u[0])
    ny = len(grid_data_u)
    X, Y = np.meshgrid(np.arange(0, 1, 1.0 / nx), np.arange(0, float(ny)/nx, 1.0 / nx))
    U = grid_data_u
    V = grid_data_v
    fig, ax = plt.subplots()
    if not has_xtick:
        ax.set_xticks(())
        ax.set_xticklabels(())
    if not has_ytick:
        ax.set_yticks(())
        ax.set_yticklabels(())
    ax.set_aspect('equal')
    plt.quiver(X, Y, U, V)
    plt.savefig(output_filename)
    plt.close(fig)
    print ('Rendered <%s>' % output_filename)