def main():
    # Create the grid
    x = arange(-100, 101)
    y = arange(-100, 101)

    # Create the meshgrid
    Y, X = meshgrid(x, y)
    A = 1
    B = 2
    V = 6*pi / 201
    W = 4*pi / 201
    F = A*sin(V*X) + B*cos(W*Y)
    Fx = V*A*cos(V*X)
    Fy = W*B*-sin(W*Y)

    # Show the images
    show_image(F)
    show_image(Fx)
    show_image(Fy)

    # Create the grid for the quivers
    xs = arange(-100, 101, 10)
    ys = arange(-100, 101, 10)

    # Here we determine the direction of the quivers
    Ys, Xs = meshgrid(ys, xs)
    FFx = V*A*cos(V*Xs)
    FFy = W*B*-sin(W*Ys)

    # Draw the quivers and the image
    clf()
    imshow(F, cmap=cm.gray, extent=(-100, 100, -100, 100))
    quiver(ys, xs, -FFy, FFx, color='red')
    show()

def computeKOPgrid():
    #compute KOP grid
    KOPgrid = pb.zeros((len(noise2),len(noise3),len(x0s)))
    i,j,k=0,0,0
    for i in range(len(noise2)):
        for j in range(len(noise3)):
            for k in range(len(x0s)):
                radical="_%i_HMR_%i_ML_CpES1_100_CpES2_100_x0_%i_noise3_%i_noise2_%i_noise1_%i_15s" %(nbn1, nbn2, int(x0s[k]*10), int(noise3[j]), int(noise2[i]*100),  int(noise2[i]*200),)
                KOP1 = np.load('traces_may22/KOP1'+radical+".npy")
                KOP2 = np.load('traces_may22/KOP2'+radical+".npy")
                KOPgrid[i,j,k] = KOP1[0,:].mean() + KOP2[0,:].mean()
    
    kopmin = KOPgrid.min()
    kopmax = KOPgrid.max()
    print KOPgrid
    #plot
    fig = pb.figure()
    ax = fig.add_subplot(311, projection='3d')
    #k=0 
    X,Y = pb.meshgrid(noise2,noise3)
    for k in range(len(x0s)):
        #img = ax.contourf(KOPgrid[:,:,k], vmin=kopmin, vmax=kopmax), offset=x0s[k])
        img = ax.contourf(X, Y, KOPgrid[:,:,k].T, zdir='z', levels=pb.linspace(kopmin, kopmax, 50), vmin=kopmin, vmax=kopmax, offset=x0s[k])
        #cbar = pb.colorbar(img)
        #fig.savefig("param3Dspace_CpES1_%i_CpES2_%i_x0_%i.png" %(int(CpES1s[i]*100), int(CpES2s[j]*100), int(x0s[k]*10)))
    ax.set_xlabel("noise2")
    ax.set_xticks(pb.linspace(min(noise2),max(noise2),5))
    ax.set_ylabel("noise3")
    ax.set_yticks(pb.linspace(min(noise3),max(noise3),5))
    ax.set_zlabel("x0")
    ax.set_zticks(pb.linspace(min(x0s),max(x0s),5))
    ax.set_zlim(min(x0s), max(x0s))
    
    ax2=fig.add_subplot(312)
    index_noise3=0;
    X2,Y2=pb.meshgrid(x0s, noise2)
    img2 = ax2.contourf(X2, Y2,KOPgrid[:,index_noise3,:],levels=pb.linspace(kopmin, kopmax, 50), vmin=kopmin, vmax=kopmax)
#    img = ax2.imshow(KOPgrid[index_noise2,:,:].T, interpolation="nearest", origin="lower")
    ax2.set_xlabel("|x0|")
    ax2.set_ylabel("noise_1_2")
    cbar=pb.colorbar(img2)
    
    noiseGrid=pb.zeros((len(noise2),len(x0s)))
    for n in range(len(noise2)):
        noiseGrid[n,:] = KOPgrid[n,0,:]
    ngmin=noiseGrid.min()
    ngmax=noiseGrid.max()
    
    ax3=fig.add_subplot(313)
    X3,Y3=pb.meshgrid(x0s, range(len(noise2)))
    img3 = ax3.contourf(X3, Y3, noiseGrid, levels=pb.linspace(ngmin, ngmax, 50), vmin=ngmin, vmax=ngmax)
    #img2 = ax3.imshow(noiseGrid.T, interpolation="nearest", origin="lower")
    cbar2=pb.colorbar(img3)
    ax3.set_xlabel("|x0|")
    ax3.set_ylabel("noise_1_2_3")
    
    print noiseGrid.T
    
    pb.show()

    def boundary_layer(self): 
        self.dpi3 = 102
        self.fig2 = Figure((3.0, 3.0), dpi=self.dpi3)
        self.axes3 = self.fig2.add_subplot(111)
        self.axes3.set_axis_bgcolor('white')
        self.x = arange(0,10,.5)
        self.y = arange(0,10,.5) 
        self.X, self.Y = pylab.meshgrid(self.x,self.y)
        self.Z, self.Zt = pylab.meshgrid(self.x,self.y)

        if not self.update:
            self.axes3.contourf(self.X,self.Y,self.Z)

    def __init__( self, t, x, y, A=[1., 0.85], a=[0.25, 0.85] ):
        '''
        
        Initializing generalized thalamo-cortical loop. Full
        functionality is only obtained in the subclasses, like DOG,
        full_eDOG, etc.
        
        Parameters
        ----------
        
        t : array
            1D Time vector
        x : np.array
            1D vector for x-axis sampling points
        y : np.array
            1D vector for y-axis sampling points

        Keyword arguments
        -----------------
        
        A : sequence (default A = [1., 0.85])
            Amplitudes for DOG receptive field for center and surround, 
            respectively
        a : sequence (default a = [0.25, 0.85])
            Width parameter for DOG receptive field for center and surround,
            respectively

        Usage
        -----
            Look at subclasses for example usage            
        '''
        # Set parameteres as attributes
        self.name = 'pyDOG Toolbox'
        self.t = t
        self.A = A
        self.a = a
        self.x = x
        self.y = y
        # Find sampling rates and sampling freqs and 
        self.nu_xs = 1./(x[1]-x[0])
        self.nu_ys = 1./(y[1]-y[0])
        self.fs = 1./(t[1]-t[0])
        self.f = fft.fftfreq(pl.asarray(t).size, t[1]-t[0])
        # fftshift spatial frequency,
        self.nu_x = fft.fftfreq(pl.asarray(x).size, x[1]-x[0])
        self.nu_y = fft.fftfreq(pl.asarray(y).size, y[1]-y[0])
        # Make meshgrids, may come in handy
        self._xx, self._yy = pl.meshgrid(self.x, self.y)
        self._nu_xx, self._nu_yy = pl.meshgrid(self.nu_x, self.nu_y)
        # r is needed for all circular rfs
        self.r = pl.sqrt(self._xx**2 + self._yy**2)
        self.k = 2 * pl.pi * pl.sqrt(self._nu_xx**2 + self._nu_yy**2)

    def init_contour(self):
        self.data2 = []
        self.x = arange(0,10,.5)
        self.y = arange(0,10,.5) 
        self.X, self.Y = pylab.meshgrid(self.x,self.y)
        self.Z, self.Zt = pylab.meshgrid(self.x,self.y)
        self.dpi2 = 101
        self.fig1 = Figure((3.0, 3.0), dpi=self.dpi2)
        self.axes2 = self.fig1.add_subplot(111)
        self.axes2.set_axis_bgcolor('white')

        if not self.update:
            initplot = self.plot_data2 = self.axes2.contourf(self.X,self.Y,self.Z)
            self.cb = self.fig1.colorbar(initplot)

    def vecPlot(self, row, col, mWx, mWy, frmSize, winSize):
        pl.clf()
        pl.ion()

        #pWy = mWy[row][col]
        #まずはx方向のWx
        pWx = mWx[row,col:(frmSize-winSize+1)+row,0,:]
        #print "pWx sum:",np.sum(pWx[0,:])
        print np.mean(pWx*pWx, axis=0)
        pWx = np.sqrt(pWx * pWx)
        r,c = pWx.shape
        x = pl.arange(c+1)
        y = pl.arange(r+1)
        X, Y = pl.meshgrid(x, y)

        pl.subplot2grid((2,2),(0,0))
        pl.pcolor(X, Y, pWx)
        pl.xlim(0,c)
        pl.ylim(0,r)
        pl.colorbar()
        pl.title("user_1 (t:"+str(row)+")")
        pl.gray()

        #いまだけ
        pWy = mWy[row,col:(frmSize-winSize+1)+row,0,:]
        #print "pWy sum:",np.sum(pWy[0,:],axis=1)
        #print np.sum(pWy*pWy,axis=1)
        #print np.mean(pWy*pWy, axis=0)
        pWy = np.sqrt(pWy * pWy)
        r,c = pWx.shape
        x = pl.arange(c+1)
        y = pl.arange(r+1)
        X, Y = pl.meshgrid(x, y)

        pl.subplot2grid((2,2),(0,1))
        pl.pcolor(X, Y, pWy)    
        pl.xlim(0,c)
        pl.ylim(0,r)
        pl.colorbar()
        pl.title("user_2 (t:"+str(col)+")")
        pl.gray()

        pl.subplot2grid((2,2),(1,0),colspan=2)
        pl.plot(np.mean(pWx*pWx, axis=0),color="r")
        pl.plot(np.mean(pWy*pWy, axis=0),color="b")

        pl.draw()

        print "pWx shape:",pWx.shape

 def __init__( self, t, x, y, contrast=1 ):
     self.t = t
     self.x = x
     self.y = y
     self.contrast = float(contrast)
     self._xx, self._yy = pl.meshgrid(self.x, self.y)
     self.make_frames()

def test_operation_approx():
    def flux_qubit_potential(phi_m, phi_p):
        return 2 + alpha - 2 * pl.cos(phi_p)*pl.cos(phi_m) - alpha * pl.cos(phi_ext - 2*phi_p)
    alpha = 0.7
    phi_ext = 2 * np.pi * 0.5
    phi_m = pl.linspace(0, 2*np.pi, 100)
    phi_p = pl.linspace(0, 2*np.pi, 100)
    X,Y = pl.meshgrid(phi_p, phi_m)
    Z = flux_qubit_potential(X, Y).T

    # the diagram creatinos
    from diagram.operations.computations import multiply
    from diagram.ternary import AEV3DD

    aevdd = AEV3DD()
    diagram3 = aevdd.create(Z, 0, True)
    diagram4 = aevdd.create(Z, 0, True)

    aevdd_mat = multiply(diagram3, diagram4, 9).to_matrix(77, True)
    aevdd_mat_approx = multiply(diagram3, diagram4, 9, approximation_precision=1, in_place='1').to_matrix(27, True)

    pl.plt.figure()
    fig, ax = pl.plt.subplots()
    p = ax.pcolor(X/(2*pl.pi), Y/(2*pl.pi), Z, cmap=pl.cm.RdBu, vmin=abs(Z).min(), vmax=abs(Z).max())
    cb = fig.colorbar(p, ax=ax)
    p = ax.pcolor(X/(2*pl.pi), Y/(2*pl.pi), Z, cmap=pl.cm.RdBu, vmin=abs(Z).min(), vmax=abs(Z).max())
    cb = fig.colorbar(p, ax=ax)
    # cnt = ax.contour(Z, cmap=pl.cm.RdBu, vmin=abs(Z).min(), vmax=abs(Z).max(), extent=[0, 1, 0, 1])
    pl.show()

def plot_thresholds(rawdata, scan_values, plane='horizontal',
                    xlabel='turns', ylabel='intensity [particles]', zlabel='normalized emittance',
                    xlimits=((0.,8192)), ylimits=((0.,7.1e11)), zlimits=((0., 10.))):

    # Prepare input data.
    # x axis
    t = rawdata[0,:,:]
    turns = plt.ones(t.shape).T * plt.arange(len(t))
    turns = turns.T

    # z axis
    epsn_abs = {}
    epsn_abs['horizontal'] = plt.absolute(rawdata[11,:,:])
    epsn_abs['vertical']   = plt.absolute(rawdata[12,:,:])

    # Prepare plot environment.
    ax11, ax13 = _create_axes(xlabel, ylabel, zlabel, xlimits, ylimits, zlimits)
    cmap = plt.cm.get_cmap('jet', 2)
    ax11.patch.set_facecolor(cmap(range(2))[-1])
    cmap = plt.cm.get_cmap('jet')

    x, y = plt.meshgrid(turns[:,0], scan_values)
    z = epsn_abs[plane]

    threshold_plot = ax11.contourf(x, y, z.T, levels=plt.linspace(zlimits[0], zlimits[1], 201),
                                   vmin=zlimits[0], vmax=zlimits[1], cmap=cmap)
    cb = plt.colorbar(threshold_plot, ax13, orientation='vertical')
    cb.set_label(zlabel)

    plt.tight_layout()

def dist2(x):
  R, GSIGMAS = pylab.meshgrid(r[r<fit_rcutoff], gsigmas)
  g = pylab.zeros_like(GSIGMAS)
  g = evalg(x, GSIGMAS, R)
  gfit = pylab.reshape(g, len(eta)*len(r[r<fit_rcutoff]))
  return gfit - pylab.reshape([g[r<fit_rcutoff] for g in ghs],
                              len(gsigmas)*len(r[r<fit_rcutoff]))

def dist2(x):
  R, ETA = pylab.meshgrid(r[r<fit_rcutoff], eta)
  g = pylab.zeros_like(ETA)
  g = evalg(x, ETA, R)
  gfit = pylab.reshape(g, len(eta)*len(r[r<fit_rcutoff]))
  return gfit - pylab.reshape([g[r<fit_rcutoff] for g in ghs],
                              len(eta)*len(r[r<fit_rcutoff]))

    def fresnelConvolutionTransform(self,d) :
        # make intensity distribution
        i2 = Intensity2D(self.nx,self.startx,self.endx,
                         self.ny,self.starty,self.endy,
                         self.wl)       

        # FT on inital distribution 
        u1ft = pl.fft2(self.i)

        # 2d convolution kernel
        k = 2*pl.pi/i2.wl
        
        # make spatial frequency matrix
        maxsfx = 2*pl.pi/self.dx
        maxsfy = 2*pl.pi/self.dy
        
        dsfx = 2*maxsfx/(self.nx)
        dsfy = 2*maxsfy/(self.ny)
        
        self.sfx = pl.arange(-maxsfx/2,maxsfx/2+1e-15,dsfx/2)
        self.sfy = pl.arange(-maxsfy/2,maxsfy/2+1e-15,dsfy/2)

        [self.sfxgrid, self.sfygrid] = pl.fftshift(pl.meshgrid(self.sfx,self.sfy))
                
        # make convolution kernel 
        kern = pl.exp(1j*d*(self.sfxgrid**2+self.sfygrid**2)/(2*k))
        
        # apply convolution kernel and invert
        i2.i = pl.ifft2(kern*u1ft) 

        return i2

def main():

    nx = floor(roh[0]/(box[1]/N))
    ny = floor(roh[1]/(box[1]/N))
   # print nx, ny
    field=zeros([N,N])
    #field[nx][ny] = 1./(dx*dx)
    field[nx+9][ny+9] = 1./(dx*dx)
    field[nx-9][ny-9] = 1./(dx*dx)
    field[nx-9][ny+9] = 1./(dx*dx)
    field[nx+9][ny-9] = 1./(dx*dx)
    #plot the initial field;
    initFig=p.figure()
    gx,gy=p.meshgrid(gridpnts,gridpnts)
    ax = p3.Axes3D(initFig)
    ax.scatter3D(ravel(gy),ravel(gx),ravel(field))
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_zlabel('$\psi(x,y)$')
    plt.savefig("ROHxySP3.png")
    
    
    
    jacobi(field)
    gauss_seidel(field)

    def plot_wire_surface_pcolor(self):
        """
        Plot the fraction of executions as a function of the fraction of nodes for
        each source. Three plots are created: wireframe, surface, and pseudocolor.

        """
        logging.debug('')
        if self.dimension > 0:
            return
        array = self._data_to_array()

        x = pylab.linspace(0.0, 1.0, len(array[0])+1)
        y = pylab.linspace(0.0, 1.0, len(array)+1)
        X, Y = pylab.meshgrid(x, y)

        #fig_wire = MyFig(self.options, xlabel='Probability p', ylabel='Fraction of Nodes', ThreeD=True)
        #fig_surf = MyFig(self.options, xlabel='Probability p', ylabel='Fraction of Nodes', ThreeD=True)
        fig_pcol = MyFig(self.options, xlabel='Probability p', ylabel='Fraction of Nodes')

        #fig_wire.ax.plot_wireframe(X, Y, array)
        #fig_surf.ax.plot_surface(X, Y, array, rstride=1, cstride=1, linewidth=1, antialiased=True)
        pcolor = fig_pcol.ax.pcolor(X, Y, array, cmap=cm.jet, vmin=0.0, vmax=1.0)
        cbar = fig_pcol.fig.colorbar(pcolor, shrink=0.8, aspect=10)
        cbar.ax.set_yticklabels(pylab.linspace(0.0, 1.0, 11), fontsize=0.8*self.options['fontsize'])

        #for ax in [fig_wire.ax, fig_surf.ax]:
            #ax.set_zlim3d(0.0, 1.01)
            #ax.set_xlim3d(0.0, 1.01)
            #ax.set_ylim3d(0.0, 1.01)

        #self.figures['wireframe'] = fig_wire.save('wireframe_' + str(self.data_filter))
        #self.figures['surface'] = fig_surf.save('surface_' + str(self.data_filter))
        self.figures['pcolor'] = fig_pcol.save('pcolor_' + str(self.data_filter))

def griddata( X, Y, Z, xl, yl, xr, yr, dx):
    # define grid.
    xi, yi = p.meshgrid( p.linspace(xl,xr, int((xr-xl)/dx)+1), p.linspace(yl,yr, int((yr-yl)/dx)+1))
    # grid the data.
    zi = mgriddata(X,Y,Z,xi,yi)
    New = grid( zi, xl, yl, dx)
    return New

    def prepare_sensors(sensor_locations, resolution=51):
        """
        Common method, to pre-process sensors before display (project them in 2D).
        """


        def sphere_fit(params):
            """Function to fit the sensor locations to a sphere"""
            return ((sensor_locations[:, 0] - params[1]) ** 2 + (sensor_locations[:, 1] - params[2]) ** 2
                    + (sensor_locations[:, 2] - params[3]) ** 2 - params[0] ** 2)


        (radius, circle_x, circle_y, circle_z) = leastsq(sphere_fit, (1, 0, 0, 0))[0]
        # size of each square
        ssh = float(radius) / resolution         # half-size
        # Generate a grid and interpolate using the gridData module
        x_arr = numpy.arange(circle_x - radius, circle_x + radius, ssh * 2.0) + ssh
        y_arr = numpy.arange(circle_y - radius, circle_y + radius, ssh * 2.0) + ssh
        x_arr, y_arr = pylab.meshgrid(x_arr, y_arr)

        # project the sensor locations onto the sphere
        sproj = sensor_locations - numpy.array((circle_x, circle_y, circle_z))
        sproj = radius * sproj / numpy.c_[numpy.sqrt(numpy.sum(sproj ** 2, axis=1))]
        sproj += numpy.array((circle_x, circle_y, circle_z))
        return dict(sproj=sproj, x_arr=x_arr, y_arr=y_arr,
                    circle_x=circle_x, circle_y=circle_y, rad=radius)

def plot_layered_images(fig):
    def func3(x,y):
        return (1- x/2 + x**5 + y**3)*exp(-x**2-y**2)

    # make these smaller to increase the resolution
    dx, dy = 0.05, 0.05

    x = arange(-3.0, 3.0, dx)
    y = arange(-3.0, 3.0, dy)
    X,Y = meshgrid(x, y)

    # when layering multiple images, the images need to have the same
    # extent.  This does not mean they need to have the same shape, but
    # they both need to render to the same coordinate system determined by
    # xmin, xmax, ymin, ymax

    xmin, xmax, ymin, ymax = min(x), max(x), min(y), max(y)
    extent = xmin, xmax, ymin, ymax
    Z1 = array(([0,1]*4 + [1,0]*4)*4); Z1.shape = 8,8  # chessboard
    Z2 = func3(X, Y)

    axes = fig.gca()
    axes.imshow(Z1, cmap=cm.gray, interpolation='nearest',
                 extent=extent)
    axes.hold(True)
    axes.imshow(Z2, cmap=cm.jet, alpha=.9, interpolation='bilinear',
                 extent=extent)

    def updateColorTable(self, cItem):
        print "now viz!"+str(cItem.row())+","+str(cItem.column())

        row = cItem.row()
        col = cItem.column()

        pl.clf()
        #pl.ion()
        x = pl.arange(self.dataDimen+1)
        y = pl.arange(self.dataDimen+1)
        X, Y = pl.meshgrid(x, y)
        pl.subplot(1,2,1)
        pl.pcolor(X, Y, self.mWx[row*self.dataMaxRange+col])
        pl.gca().set_aspect('equal')
        pl.colorbar()
        pl.gray()
        pl.title("user 1")

        pl.subplot(1,2,2)
        pl.pcolor(X, Y, self.mWy[row*self.dataMaxRange+col])
        pl.gca().set_aspect('equal')
        pl.colorbar()
        pl.gray()
        pl.title("user 2")
        #pl.tight_layout()

        pl.draw()
        #pl.show()
        pl.show(block=False) 

def demo3():
    theta = RiemannTheta_Function()
    Omega = np.array([[1.j, .5], [.5, 1.j]])
    print
    print "Calculating 1,000,000 points of the Riemann Theta Function on GPU..."
    print
    print "Omega = [i   .5]"
    print "        [.5   i]"
    print 
    print "z = (x + iy, 0) where (0 < x < 1) and (0 < y < 5)"
    SIZE = 1000
    x = np.linspace(0,1,SIZE)
    y = np.linspace(0,5,SIZE)
    X,Y = p.meshgrid(x,y)
    Z = X + Y*1.0j
    Z = Z.flatten()
    Z = [[z,0] for z in Z]
    print "Starting computation on the GPU"
    start = time.clock()
    U,V = theta.exp_and_osc_at_point(Z, Omega, batch=True)
    done = time.clock() - start
    print "Time to perform the calculation: " + str(done)
    print
    print "Starting computation on the CPU"
    start = time.clock()
    U,V = theta.exp_and_osc_at_point(Z, Omega, batch=True, gpu=False)
    done = time.clock() - start
    print "Time to perform the calculation: " + str(done)
    print

def plotInit(Plotting, Elements):
	if (Plotting == 2):
		loc = [i.xy for i in Elements]
		x = [i.real for i in loc]
		y = [i.imag for i in loc]
		x = list(sorted(set(x))) 
		x.remove(-10)
		y = list(sorted(set(y)))

		X, Y = pylab.meshgrid(x, y)
		U = pylab.ones(shape(X))
		V = pylab.ones(shape(Y))

		pylab.ion()
		fig, ax = pylab.subplots(1,1)
		graph = ax.quiver(X, Y, U, V)
		pylab.draw()
	else:
		pylab.ion()
		graph, = pylab.plot(1, 'ro', markersize = 2) 
		x = 2
		pylab.axis([-x,x,x,-x])

		graph.set_xdata(0)
		graph.set_ydata(0)
		pylab.draw()

	return graph