def baltimore_gmm(data):
    def fgmm(x):
        return abs(np.sum(gmmimmat[gmmimmat > x]) * .0005 ** 2 - 0.95)

    model = gmm.GMM(3)
    model.train(data, random=False)

    X, Y = makegrid(data)
    gmmimmat = np.zeros(X.shape)

    for i in xrange(X.shape[0]):
        for j in xrange(X.shape[1]):
            gmmimmat[i, j] = model.dgmm(np.array([X[i, j], Y[i, j]]))

    plt.jet()
    plt.imshow(gmmimmat, origin='lower')
    plt.ylim([0, X.shape[0]])
    plt.xlim([0, X.shape[1]])
    plt.savefig('baltimore_gmm.pdf')

    thresh = opt.fmin(fgmm, 10)[0]
    bools = gmmimmat > thresh
    mat = np.zeros(X.shape)
    mat += bools
    plt.imshow(mat, origin='lower')
    plt.ylim([0, X.shape[0]])
    plt.xlim([0, X.shape[1]])
    plt.savefig('gmmavoid.pdf')

def plotear(xi,yi,zi):
    # mask inner circle
    interior1 = sqrt(((xi+1.5)**2) + (yi**2)) < 1.0 
    interior2 = sqrt(((xi-1.5)**2) + (yi**2)) < 1.0
    zi[interior1] = ma.masked
    zi[interior2] = ma.masked
    p.figure(figsize=(16,10))
    pyplot.jet()
    max=2.8
    min=0.4
    steps = 50
    levels=list()
    labels=list()
    for i in range(0,steps):
	levels.append(int((max-min)/steps*100*i)*0.01+min)
    for i in range(0,steps/2):
	labels.append(levels[2*i])
    CSF = p.contourf(xi,yi,zi,levels,norm=colors.LogNorm())
    CS = p.contour(xi,yi,zi,levels, format='%.3f', labelsize='18')
    p.clabel(CS,labels,inline=1,fontsize=9)
    p.title('electrostatic potential of two spherical colloids, R=lambda/3',fontsize=24)
    p.xlabel('z-coordinate (3*lambda)',fontsize=18)
    p.ylabel('radial coordinate r (3*lambda)',fontsize=18)
    # add a vertical bar with the color values
    cbar = p.colorbar(CSF,ticks=labels,format='%.3f')
    cbar.ax.set_ylabel('potential (reduced units)',fontsize=18)
    cbar.add_lines(CS)
    p.show()

def plot_percent_traffic(category, num_bins, city):
    '''
    plots percent traffic (review_counts) of category in bin
    '''
    a, b =bin_by_review_count(category, num_bins, city)
    c = []
    max_ratio = 0
    for i in xrange(len(a)):
        c.append([])
        for j in xrange(len(a[i])):
            try: 
                ratio = a[i][j]/float(b[i][j])
                c[i].append(ratio)
                if ratio > max_ratio: 
                    max_ratio = ratio
#                     print i,j,max_ratio
            except ZeroDivisionError:
                c[i].append('NA')
#     print c
    
    for i in xrange(len(c)):
        for j in xrange(len(c[i])):
            if c[i][j] == 'NA':
                c[i][j] = 0
                
    plt.imshow(c, interpolation='none', alpha = 1)

    max_lat, min_lat, max_lon, min_lon = city_edges(city)
    plt.xticks([0,len(c[0])-1],[min_lon, max_lon])
    plt.yticks([0,len(c)-1],[min_lat, max_lat])

    plt.jet()
    cb = plt.colorbar() #make color bar
    cb.set_ticks([0,max_ratio])   #two ticks
    cb.set_ticklabels(['low traffic','lots of traffic'])  # put text labels on them

def draw(ax, mu, Sigma):
    # 描画のクリア
    ax.collections = []

    # 等高線を描画
    xlist = np.linspace(-5, 5, 50)
    ylist = np.linspace(-5, 5, 50)
    x,y = np.meshgrid(xlist,ylist)
    z = np.zeros((50,50))
    for i in range(len(ylist)):
        for j in range(len(xlist)):
            xx = np.array([[xlist[j]], [ylist[i]]])
            z[i,j] = gaussian(xx, mu, Sigma)
    cs = ax.pcolor(x, y, z)
    plt.colorbar(cs)
    plt.jet()
    #plt.bone()

    ax.contour(x, y, z, np.linspace(0.0001,0.5,25), colors='k', linewidth=1)

    # ガウス分布の平均を描画
    ax.scatter(mu[0], mu[1], c='b', marker='x')

    # 軸の調整
    ax.set_xlim(-5,5)
    ax.set_ylim(-5,5)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title('2D Normal distribution')
    ax.grid()

    plt.draw()

def plot_bins(category,num_bins, city):
    '''
    takes a category and city and plots the bins for it
    '''
    a, b = make_bins(category, num_bins, city)
    c = []
    max_ratio = 0
    for i in xrange(len(a)):
        c.append([])
        for j in xrange(len(a[i])):
            try: 
                ratio = a[i][j]/float(b[i][j])
                c[i].append(ratio)
                if ratio > max_ratio: 
                    max_ratio = ratio
#                     print i,j,max_ratio
            except ZeroDivisionError:
                c[i].append('NA')
#     print c
    
    for i in xrange(len(c)):
        for j in xrange(len(c[i])):
            if c[i][j] == 'NA':
                c[i][j] = max_ratio
                
    plt.imshow(c, interpolation='none', alpha = 1)

    max_lat, min_lat, max_lon, min_lon = city_edges(city)
    plt.xticks([0,len(c[0])-1],[min_lon, max_lon])
    plt.yticks([0,len(c)-1],[min_lat, max_lat])

    plt.jet()
    cb = plt.colorbar() #make color bar
    cb.set_ticks([max_ratio, 0])   #two ticks
    cb.set_ticklabels(['high concentration', 'low concentration'])  # put text labels on them

def plot_saccade_stats(sac,bins=36, fig=None):
    """
    Draws individual saccades (polar coordinatess), directional histogram,
    individual saccades (cartesian coordinates), and a histogram of saccade
    peak velocities
    """
    if fig is None:
        fig = plt.figure()
    dx = sac.xf - sac.xi
    dy = sac.yf - sac.yi
    radii = sac.amplitude
    thetas = np.arctan2(dy, dx)
    theta_bins = np.linspace(-np.pi,np.pi,bins+1)
    theta_hist = np.histogram(thetas,bins=theta_bins)

    plt.jet()
    plt.title("individual saccades")
    # XXX: there's an interesting artifact if we use sac.amplitude as the
    # sac_length in the scatter plot below
    sac_length = np.sqrt((dx**2+dy**2))
    plt.subplot(221,polar=True)
    plt.scatter(thetas,sac_length,alpha=.5,c=sac.amplitude,s=sac.vpeak/10.0)
    plt.colorbar()
    plt.subplot(222,polar=True)
    plt.title("Directional histogram")
    bars = plt.bar(theta_bins[:-1],theta_hist[0],width=(2*np.pi)/bins, alpha=.5)
    plt.subplot(223)
    plt.scatter(dx,dy,alpha=.5,c=sac.amplitude,s=sac.vpeak/10.0)

    #plt.hist(sac_length,bins=100)
    #plt.xlabel("saccade lengths (pixels)")
    plt.subplot(224)
    plt.hist(sac.vpeak,bins=100)
    plt.xlabel("saccade peak velocities")

def test_pca_white(sh=(12, 12), m=500, eps=0.1, rc=(10, 10), out_dir="img"):
    # Visualize examples in r x c grid of pca whitened images.
    # Also show that the covariance matrix of whitening data
    # is a diagonal matrix with descending values
    x = test_pca_input(sh, m)
    pca = PCA(x)

    x_white, mean = pca.whiten_pca(x, eps=eps)

    rows, cols = rc

    f, axes = plt.subplots(rows, cols, sharex='col', sharey='row')

    plt.subplots_adjust(hspace=0.1, wspace=0)
    plt.jet()

    for r in range(rows):
        for c in range(cols):
            axes[r][c].imshow(x_white[:, r*cols + c].reshape(sh))
            axes[r][c].axis('off')

    path = os.path.join(out_dir, "pca_filters_{}_{}_{}.png".format(sh, m, eps))
    print("saving {}...".format(path))
    plt.savefig(path, bbox_inches='tight')

    cv = np.dot(x_white, x_white.T)
    plt.clf()
    plt.imshow(cv)
    plt.show()

def image_gen(dataName, initValue, finalValue, increment,imgdpi):
	for i in range(initValue,finalValue,increment):
		if not os.path.exists(dataName+"r_"+str(i)+"_abspsi2.png"):
			real=open(dataName + '_' + str(i)).read().splitlines()
			img=open(dataName + 'i_' + str(i)).read().splitlines()
			a_r = numpy.asanyarray(real,dtype='f8') #64-bit double
			a_i = numpy.asanyarray(img,dtype='f8') #64-bit double
			a = a_r[:] + 1j*a_i[:]
			b = np.reshape(a,(xDim,yDim))
			f = plt.imshow(abs(b)**2)
			plt.jet()
			plt.gca().invert_yaxis()
			plt.savefig(dataName+"r_"+str(i)+"_abspsi2.png",dpi=imgdpi)
			plt.close()
			g = plt.imshow(np.angle(b))
			plt.gca().invert_yaxis()
			plt.savefig(dataName+"r_"+str(i)+"_phi.png",dpi=imgdpi)
			plt.close()
			f = plt.imshow(abs(np.fft.fftshift(np.fft.fft2(b)))**2)
			plt.gca().invert_yaxis()
			plt.jet()
			plt.savefig(dataName+"p_"+str(i)+"_abspsi2.png",dpi=imgdpi)
			plt.close()
			g = plt.imshow(np.angle(np.fft.fftshift(np.fft.fft2(b))))
			plt.gca().invert_yaxis()
			plt.savefig(dataName+"p_"+str(i)+"_phi.png",dpi=imgdpi)
			plt.close()
			print "Saved figure: " + str(i) + ".png"
			plt.close()
		else:
			print "File(s) " + str(i) +".png already exist."

def makeConfMat(estClasses, gtClasses, outFilename, numClasses = None, plotLabels = False):
   #If not defined, find number of unique numbers in gtClasses
   if numClasses == None:
      numClasses = len(np.unique(gtClasses))

   #X axis is est, y axis is gt
   #First index is y, second is x
   confMat = np.zeros((numClasses, numClasses))
   numInstances = len(estClasses)
   for (gtIdx, estIdx) in zip(gtClasses.astype(int), estClasses.astype(int)):
      confMat[gtIdx, estIdx]  += 1

   plt.jet()
   plt.matshow(confMat)
   plt.colorbar()
   plt.xlabel("Est class")
   plt.ylabel("True class")
   plt.title("Confusion matrix")
   ax = plt.gca()
   ax.xaxis.set_ticks_position('bottom')

   #Plot labels for each field
   if plotLabels:
      for i in range(numClasses):
         for j in range(numClasses):
            labelStr = generateStatString(confMat, i, j)
            #text receives x, y coord of plot
            ax.text(j, i, labelStr, fontweight='bold',
                  horizontalalignment='center', verticalalignment='center',
                  bbox={'facecolor':'white'}, fontsize=6)

   #plt.show()
   plt.savefig(outFilename)

def baltimore_kde(data):
    def fkde(x):
        return abs(np.sum(kdeimmat[kdeimmat > x]) * .0005 ** 2 - 0.95)

    X, Y = makegrid(data)
    kdeimmat = np.zeros(X.shape)
    kernel = stats.gaussian_kde(data.T)
    for i in xrange(X.shape[0]):
        for j in xrange(X.shape[1]):
            kdeimmat[i, j] = kernel.evaluate(np.array([X[i, j], Y[i, j]]))

    plt.jet()
    plt.imshow(kdeimmat, origin='lower')
    plt.ylim([0, X.shape[0]])
    plt.xlim([0, X.shape[1]])
    plt.savefig('baltimore_kde.pdf')

    thresh = opt.fmin(fkde, 10)[0]
    bools = kdeimmat > thresh
    mat = np.zeros(X.shape)
    mat += bools
    plt.imshow(mat, origin='lower')
    plt.ylim([0, X.shape[0]])
    plt.xlim([0, X.shape[1]])
    plt.savefig('kdeavoid.pdf')

def struct_fact(density,name,imgdpi):
	fig, ax = plt.subplots()
	#f = plt.quiver(gx,gy)
	f = plt.imshow((np.abs(np.fft.fftshift(np.fft.fft2(density)))),cmap=plt.get_cmap('prism'))
	cbar = fig.colorbar(f)
	cbar.set_clim(1e6,1e11)
	plt.jet()
	plt.savefig(name + "_struct_log10.png",dpi=imgdpi)
	plt.close()

def heatmap_plot(data, size, ratio, dir_name):
	im = plt.imshow(data, interpolation='none', aspect=ratio) # change the aspect if needed
	plt.xticks(range(size))
	plt.jet()
	plt.colorbar()
	plt.clim(0,BAR_RANGE)
	# plt.show()
	plt.savefig(dir_name + 'comp.png')
	plt.close()

def scaleAxis(data,dataName,label,value,imgdpi):
	fig, ax = plt.subplots()
	ax.xaxis.set_major_locator(ScaledLocator(dx=dx))
	ax.xaxis.set_major_formatter(ScaledLocator(dx=dx))
	f = plt.imshow(abs(data)**2)
	cbar = fig.colorbar(f)
	plt.gca().invert_yaxis()
	plt.jet()
	plt.savefig(dataName+"r_"+str(value)+"_"+label +".png",dpi=imgdpi)
	plt.close()

def run3Dheatmap(X):
    
    plt.imshow(data, interpolation='none', aspect=3./20)

    plt.xticks(range(3), ['a', 'b', 'c'])

    plt.jet()
    plt.colorbar()

    plt.show() 

def opPot(dataName,imgdpi):
	data = open(dataName).read().splitlines()
	a = numpy.asanyarray(data,dtype='f8')
	b = np.reshape(a,(xDim,yDim))
	fig, ax = plt.subplots()
	f = plt.imshow((b))
	plt.gca().invert_yaxis()
	cbar = fig.colorbar(f)
	plt.jet()
	plt.savefig(dataName + ".png",dpi=imgdpi)
	plt.close()

def heatmap_plot(data, size_x, size_y, dir_name):
	im = plt.imshow(data, interpolation='none', aspect=1.0) # 1.0 = square cell
	plt.xlim([0.5, size_x + 0.5])
	plt.ylim([0.5, size_y + 0.5])
	plt.xticks(range(1, size_x + 1), fontsize=X_SIZE)
	plt.yticks(list(reversed(range(1, size_y + 1))), fontsize=Y_SIZE)
	plt.jet()
	cb = plt.colorbar()
	cb.ax.tick_params(labelsize=CB_SIZE) 
	plt.clim(0,BAR_RANGE)
	plt.set_cmap('gray_r')
	plt.savefig(dir_name + 'comp.png')
	plt.close()

    def plotdataarray(self):
        """Plot the image array

        return axes
        """
        self.ob_plot = self.slotfig.add_axes([0.10,0.10,0.8,0.80], autoscale_on=True)
        plt.setp(plt.gca(),xticks=[],yticks=[])
        plt.jet()
        self.array=self.struct[int(self.pid[self.id])].data
        self.imarr=self.ob_plot.imshow(self.array,origin='lower')
        #Plot the apertures
        self.cbox,=self.plotbox('#00FF00',self.cx[self.id],self.cy[self.id],self.radius,self.npoint,self.naxis1, self.naxis2)
        self.tbox,=self.plotbox('#FFFF00',self.tx[self.id],self.ty[self.id],self.radius,self.npoint,self.naxis1, self.naxis2)

def demoCDIRECT(maxiter=25):
    """
    Test and visualize cDIRECT on a 2D function.  This will draw the contours
    of the target function, the final set of rectangles and mark the optimum
    with a red dot.
    """
    import matplotlib.pyplot as plt
    
    def foo(x):
        """
        Code for the Shekel function S5.  The dimensionality 
        of x is 2.  The  minimum is at [4.0, 4.0].
        """
        # return min(-.5, -sum(1./(dot(x-a, x-a)+c) for a, c in SHEKELPARAMS))
        return sin(x[0]*2)+abs(x[0]-15) + sin(x[1])+.2*abs(x[1]-6)
    # def foo(x):
    #     # several local minima, global minimia is at bottom left
    #     return 2.5 + sin((x[0]-.4)*8)+sin((x[1]+.5)*5) + .1* sum(sin(x[0]*50)) + .1* sum(sin(x[1]*50))+ x[0]*.1 - x[1] * .1

    bounds = [(1.2, 28.), (0.1, 13.)]
    optv, optx = cdirect(foo, bounds, maxiter=maxiter)
    print '***** opt val =', optv
    print '***** opt x   =', optx
    
    plt.figure(2)
    plt.clf()
    
    # plot rectangles
    c0 = [(i/100.)*(bounds[0][1]-bounds[0][0])+bounds[0][0] for i in xrange(101)]
    c1 = [(i/100.)*(bounds[1][1]-bounds[1][0])+bounds[1][0] for i in xrange(101)]
    z = array([[foo([i, j]) for i in c0] for j in c1])
    
    ax = plt.subplot(111)
    B = [array([1.2, 0.1]), array([28., 13.])]
    for line in open('finalrecs.txt').readlines():
        dat = line.strip().split(',')
        lb = array([double(x) for x in dat[0].split()])*(B[1]-B[0])+B[0]
        ub = array([double(x) for x in dat[1].split()])*(B[1]-B[0])+B[0]
        ax.add_artist(plt.Rectangle(lb, ub[0]-lb[0], ub[1]-lb[1], fc='y', ec='k', lw=1, alpha=0.25, fill=True))

    ax.plot(optx[0], optx[1], 'ro')
    cs = ax.contour(c0, c1, z, 10)
    ax.clabel(cs)
    plt.jet()
    ax.set_xlim(bounds[0])
    ax.set_ylim(bounds[1])
    ax.set_xlabel('x[0]')
    ax.set_ylabel('x[1]')
    ax.set_title('final optimum')

def test_pca(sh=(12, 12), m=10, retain=0.9):
    """
    Show plot demonstrating that the covariance matrix of x_tilde
    is a diagonal matrix with descending values
    :param sh: shape of image patches to test
    :param m: The number of samples to load
    :param retain: Percentage of variance to retain
    """
    x = test_pca_input(sh, m)
    pca = PCA(x)

    xt, mean = pca.reduce(x, retain=retain)
    cov = np.dot(xt, xt.T)

    plt.jet()
    plt.imshow(cov)
    plt.show()

def vectorfield(show=False):
    """ Data for vectorfield """
    
    R = 7.0
    J = 100.0
    
    YY, XX = np.mgrid[-10:10,-10:10]
        
    print XX
    print YY
    
    MAG = np.sqrt(XX**2 + YY**2)
    
    FIELDX_IN = np.zeros_like(XX)
    FIELDY_IN = np.zeros_like(YY)
    FIELDX_OUT = FIELDX_IN.copy()
    FIELDY_OUT = FIELDY_IN.copy()
    
    FIELDX_IN = -1*YY*(J/2.)
    FIELDY_IN = XX*(J/2.)

    FIELDX_OUT = ((R**2*J)/(2*(XX**2 + YY**2)))*(-1*YY)
    FIELDY_OUT = ((R**2*J)/(2*(XX**2 + YY**2)))*(XX)

    FIELDX = np.zeros_like(XX)
    FIELDY = np.zeros_like(YY)
    
    MAG = np.sqrt(XX**2 + YY**2)
    
    FIELDX[MAG<R] = FIELDX_IN[MAG < R]
    FIELDX[MAG>=R] = FIELDX_OUT[MAG >= R]
    
    FIELDY[MAG<R] = FIELDY_IN[MAG < R]
    FIELDY[MAG>=R] = FIELDY_OUT[MAG >= R]
    
    np.savetxt('vectorfield_x.txt', FIELDX)
    np.savetxt('vectorfield_y.txt', FIELDY)
    
    if show:
        f = plt.figure(199)
        plt.clf()
        plt.jet()
        plt.quiver(FIELDX, FIELDY, (FIELDX**2 + FIELDY**2))