def ST_beta(r,steps = 99):
	"""
	Makes a heat map of beta over ST space, for a given value of r.
	Uses the full model.
	Steps : int {99}
		Number of steps to sample S and T at
	Returns
	=======
	None
	"""

	Ss = np.linspace(-1,3,steps)
	Ts = np.linspace(0,4,steps)
	dataFull = np.zeros( (steps,steps) )
	for i,S in enumerate(Ss):
		for j,T in enumerate(Ts):
			mod = ESS_class(r,S,T)
			state = mod.getESS()
			dataFull[i,j] = mod.beta(state,r)

	#pl.figure()
	cmap = pl.get_cmap('coolwarm')
	pl.imshow(dataFull, origin = [Ts[0], Ss[0]],\
			extent = [ Ts[0],Ts[-1],Ss[0],Ss[-1] ], vmin = -1,vmax = 1, cmap = cmap)
	pl.plot( [ Ts[0], Ts[-1] ],[ 0,0 ],color='black', linewidth = 2.5 )
	pl.plot( [ 1, 1 ],[ Ss[0],Ss[-1] ],'--',color='black',  linewidth = 2.5 )
	pl.plot( [ 0, 3 ],[ 2, -1 ],'--',color='black',  linewidth = 2.5 )

def plot_Barycenter(dataset_name, feat, unfeat, repo):

    if dataset_name==MNIST:
        _, _, test=get_data(dataset_name, repo, labels=True)
        xtest1,_,_, labels,_=test
    else:
        _, _, test=get_data(dataset_name, repo, labels=False)
        xtest1,_,_ =test
        labels=np.zeros((len(xtest1),))
    # get labels
    def bary_wdl2(index): return _bary_wdl2(index, xtest1, feat, unfeat)
    
    n=xtest1.shape[-1]
    
    num_class = (int)(max(labels)+1)
    barys=[bary_wdl2(np.where(labels==i)) for i in range(num_class)]
    pl.figure(1, (num_class, 1))
    for i in range(num_class):
        pl.subplot(1,10,1+i)
        pl.imshow(barys[i][0,0,:,:],cmap='Blues',interpolation='nearest')
        pl.xticks(())
        pl.yticks(())
        if i==0:
            pl.ylabel('DWE Bary.')
        if num_class >1:
            pl.title('{}'.format(i))
    pl.tight_layout(pad=0,h_pad=-2,w_pad=-2) 
    pl.savefig("imgs/{}_dwe_bary.pdf".format(dataset_name))

    def plot_matches(self, name, show_below = True, match_maximum = None):
        """ 対応点を線で結んで画像を表示する
          入力： im1,im2（配列形式の画像）、locs1,locs2（特徴点座標）
             machescores（match()の出力）、
             show_below（対応の下に画像を表示するならTrue）"""
        im1 = self._image_1.get_array_image()
        im2 = self._image_2.get_array_image()
        self.appendimages()
        im3 = self._append_image
        if self._match_score is None:
            self.match()
        locs1 = self._image_1.get_shift_location()
        locs2 = self._image_2.get_shift_location()
        if show_below:
            im3 = numpy.vstack((im3,im3))
        pylab.figure(dpi=160)
        pylab.gray()
        pylab.imshow(im3, aspect = 'auto')

        cols1 = im1.shape[1]
        match_num = 0
        for i,m in enumerate(self._match_score):
            if m > 0 : 
                pylab.plot([locs1[i][0],locs2[m][0]+cols1], [locs1[i][1],locs2[m][1]], 'c')
                match_num = match_num + 1
            if match_maximum is not None and match_num >= match_maximum:
                break
        pylab.axis('off')
        pylab.savefig(name, dpi=160)

def makeContourPlot(scores, average, HEIGHT, WIDTH, outputId, maskId, plt_title, outputdir, barcodeId=-1, vmaxVal=100):
    pylab.bone()
    #majorFormatter = FormatStrFormatter('%.f %%')
    #ax = pylab.gca()
    #ax.xaxis.set_major_formatter(majorFormatter)
    
    pylab.figure()
    ax = pylab.gca()
    ax.set_xlabel(str(WIDTH) + ' wells')
    ax.set_ylabel(str(HEIGHT) + ' wells')
    ax.autoscale_view()
    pylab.jet()
    
    pylab.imshow(scores,vmin=0, vmax=vmaxVal, origin='lower')
    pylab.vmin = 0.0
    pylab.vmax = 100.0
    ticksVal = getTicksForMaxVal(vmaxVal)
    pylab.colorbar(format='%.0f %%',ticks=ticksVal)
    print "'%s'" % average
    if(barcodeId!=-1):
        if(barcodeId==0): maskId = "No Barcode Match,"
        else:             maskId = "Barcode Id %d," % barcodeId
    if plt_title != '': maskId = '%s\n%s' % (plt_title,maskId)
    print "Checkpoint A"
    pylab.title('%s Loading Density (Avg ~ %0.f%%)' % (maskId, average))
    pylab.axis('scaled')
    print "Checkpoint B"
    pngFn = outputdir+'/'+outputId+'_density_contour.png'
    print "Try save to", pngFn;
    pylab.savefig(pngFn, bbox_inches='tight')
    print "Plot saved to", pngFn;

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 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 window_fn_matrix(Q,N,num_remov=None,save_tag=None,lms=None):
    Q = n.matrix(Q); N = n.matrix(N)
    Ninv = uf.pseudo_inverse(N,num_remov=None) # XXX want to remove dynamically
    #print Ninv 
    info = n.dot(Q.H,n.dot(Ninv,Q))
    M = uf.pseudo_inverse(info,num_remov=num_remov)
    W = n.dot(M,info)

    if save_tag!=None:
        foo = W[0,:]
        foo = n.real(n.array(foo))
        foo.shape = (foo.shape[1]),
        print foo.shape
        p.scatter(lms[:,0],foo,c=lms[:,1],cmap=mpl.cm.PiYG,s=50)
        p.xlabel('l (color is m)')
        p.ylabel('W_0,lm')
        p.title('First Row of Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W.pdf'.format(fig_loc,save_tag))
        p.clf()

        print 'W ',W.shape
        p.imshow(n.real(W))
        p.title('Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W_im.pdf'.format(fig_loc,save_tag))
        p.clf()


    return W

    def draw(self):

        print self.edgeno

        pos = 0
        dy = 8
        edgeno = self.edgeno
        edge = self.edges[edgeno]
        edgeprev = self.edges[edgeno-1]
        p = np.round(edge["top"](1024))
        top = min(p+2*dy, 2048)
        bot = min(p-2*dy, 2048)
        self.cutout = self.flat[1][bot:top,:].copy()

        pl.figure(1)
        pl.clf()
        start = 0
        dy = 512
        for i in xrange(2048/dy):
            pl.subplot(2048/dy,1,i+1)
            pl.xlim(start, start+dy)

            if i == 0: pl.title("edge %i] %s|%s" % (edgeno,
                edgeprev["Target_Name"], edge["Target_Name"]))


            pl.subplots_adjust(left=.07,right=.99,bottom=.05,top=.95)

            pl.imshow(self.flat[1][bot:top,start:start+dy], extent=(start,
                start+dy, bot, top), cmap='Greys', vmin=2000, vmax=6000)

            pix = np.arange(start, start+dy)
            pl.plot(pix, edge["top"](pix), 'r', linewidth=1)
            pl.plot(pix, edgeprev["bottom"](pix), 'r', linewidth=1)
            pl.plot(edge["xposs_top"], edge["yposs_top"], 'o')
            pl.plot(edgeprev["xposs_bot"], edgeprev["yposs_bot"], 'o')


            hpp = edge["hpps"]
            pl.axvline(hpp[0],ymax=.5, color='blue', linewidth=5)
            pl.axvline(hpp[1],ymax=.5, color='red', linewidth=5)

            hpp = edgeprev["hpps"]
            pl.axvline(hpp[0],ymin=.5,color='blue', linewidth=5)
            pl.axvline(hpp[1],ymin=.5,color='red', linewidth=5)


            if False:
                L = top-bot
                Lx = len(edge["xposs"])
                for i in xrange(Lx):
                    xp = edge["xposs"][i]
                    frac1 = (edge["top"](xp)-bot-1)/L
                    pl.axvline(xp,ymin=frac1)

                for xp in edgeprev["xposs"]: 
                    frac2 = (edgeprev["bottom"](xp)-bot)/L
                    pl.axvline(xp,ymax=frac2)

            start += dy

def ShowDynamicalResults(indexmap,output):
    import pylab
    pylab.figure()
    for i in range(4):
        pylab.subplot('22'+str(i+1))
        pylab.imshow(indexmap[i],interpolation='nearest')
    pylab.savefig(output+'.png')

def centerofmass(array,threshold=None,useBrightest=0):
    '''
    array is a float 32 2-dimensional numpy array
    cent = numpy.array([x,y])
    '''
    newarray = None
    useBrightest = int(useBrightest)
    if(threshold!=None):
        boolarray = array>=threshold
        newarray = boolarray*array
    else:
        newarray = array + 0

    if(useBrightest>0):
        newarray = brightest(newarray,useBrightest)

    #Classic
    cent = np.zeros(2)
    totalmass = float(newarray.sum())
    Xgrid,Ygrid = np.meshgrid(np.arange(newarray.shape[1]),np.arange(newarray.shape[0]))
    if totalmass<1.0:
        print 'totalmass: ',
        print totalmass
        imshow(newarray)
        show()

    cent[1] = np.sum(Ygrid*newarray)/totalmass
    cent[0] = np.sum(Xgrid*newarray)/totalmass
    return cent

def explore_data(data, images, target):

	# try to determine the type of data...
	print "data_type belonging to key data:"
	try: 
		print np.dtype(data)
	except TypeError as err:
		print err
	
	print "It has dimension", np.shape(data)

	# plot a 3
	
	# get indices of all threes in target
	threes = np.where(target == 3)
	#assert threes is not empty
	assert(len(threes) > 0)
	# choose the first 3
	three_indx = threes[0]
	# get the image
	img = images[three_indx][0]

	#plot it
	plot.figure()
	plot.gray()
	plot.imshow(img, interpolation = "nearest")
	plot.show()
	plot.close()

def main():
    src_cv_img_1 = cv.LoadImage("data/gorskaya/images/1/g126.jpg")
    src_cv_img_gray_1 = cv.LoadImage("data/gorskaya/images/1/g126.jpg", 
        cv.CV_LOAD_IMAGE_GRAYSCALE)

    src_cv_img_2 = cv.LoadImage("data/gorskaya/images/1/g127.jpg")
    src_cv_img_gray_2 = cv.LoadImage("data/gorskaya/images/1/g127.jpg", 
        cv.CV_LOAD_IMAGE_GRAYSCALE)

    (keypoints_1, descriptors_1) = \
        cv.ExtractSURF(src_cv_img_gray_1, None, cv.CreateMemStorage(), 
            (0, 30000, 3, 1))
    (keypoints_2, descriptors_2) = \
        cv.ExtractSURF(src_cv_img_gray_2, None, cv.CreateMemStorage(), 
            (0, 30000, 3, 1))

    print("Found {0} and {1} keypoints".format(
        len(keypoints_1), len(keypoints_2)))

    src_arr_1 = array(src_cv_img_1[:, :])[:, :, ::-1]
    src_arr_2 = array(src_cv_img_2[:, :])[:, :, ::-1]

    pylab.rc('image', interpolation='nearest')

    pylab.subplot(121)
    pylab.imshow(src_arr_1)
    pylab.plot(*zip(*[k[0] for k in keypoints_1]), 
        marker='.', color='r', ls='')

    pylab.subplot(122)
    pylab.imshow(src_arr_2)
    pylab.plot(*zip(*[k[0] for k in keypoints_2]), 
        marker='.', color='r', ls='')

    pylab.show()

def ST_pi_pure(r,steps = 99):
	"""
	Makes a heat map over ST space for a given value of r, using the pure strategy model, plotting
	fitness.

	Inputs
	======
	r : float
		Value of relatedness
	steps : int {99}
		Number of points to sample S and T at.

	"""

	Ss = np.linspace(-1,2,steps)
	Ts = np.linspace(0,3,steps)
	dataFull = np.zeros( (steps,steps) )
	for i,S in enumerate(Ss):
		for j,T in enumerate(Ts):
			x = ESS_pure(r,S,T)
			dataFull[i,j] = fit_pure(x,r,S,T)/maximal_possible_fitness(S,T)

	#pl.figure()
	cmap = pl.get_cmap('Reds')
	pl.imshow(dataFull, origin = [Ts[0], Ss[0]], interpolation = 'nearest',\
			extent = [ Ts[0],Ts[-1],Ss[0],Ss[-1] ], cmap = cmap, vmin = 0, vmax = .5*(Ss[-1]+Ts[-1]))
	pl.plot( [ Ts[0], Ts[-1] ],[ 0,0 ],color='black', linewidth = 2.5 )
	pl.plot( [ 1, 1 ],[ Ss[0],Ss[-1] ],'--',color='black',  linewidth = 2.5 )
	pl.plot( [ 0, 3 ],[ 2, -1 ],'--',color='black',  linewidth = 2.5 )

def ST_xA(r,steps = 99):
	"""
	Makes a heat map for x over ST space for a given r using the full model.
	inputs
	======
	r : float
		Value of relatedness
	steps: int {99}
		Number of decrete points at which to sample S and T

	"""


	Ss = np.linspace(-1,2,steps)
	Ts = np.linspace(0,3,steps)
	dataFull = np.zeros( (steps,steps) )
	for i,S in enumerate(Ss):
		for j,T in enumerate(Ts):
			mod = ESS_class(r,S,T)
			state = mod.getESS()
			dataFull[i,j] = mod.xA(state,r)

	#pl.figure()
	cmap = pl.get_cmap('Reds')
	pl.imshow(dataFull, origin = [Ts[0], Ss[0]],\
			extent = [ Ts[0],Ts[-1],Ss[0],Ss[-1] ], vmin = 0,vmax = 1, cmap = cmap)
	pl.plot( [ Ts[0], Ts[-1] ],[ 0,0 ],color='black', linewidth = 2.5 )
	pl.plot( [ 1, 1 ],[ Ss[0],Ss[-1] ],'--',color='black',  linewidth = 2.5 )
	pl.plot( [ 0, 3 ],[ 2, -1 ],'--',color='black',  linewidth = 2.5 )

def plotMatrix(matrix, color_scheme, row_headers, col_headers, vmin, vmax, options):

    pylab.imshow(matrix,
                 cmap=color_scheme,
                 origin='lower',
                 vmax=vmax,
                 vmin=vmin,
                 interpolation='nearest')

    # offset=0: x=center,y=center
    # offset=0.5: y=top/x=right
    offset = 0.0

    if options.xticks:
        pylab.xticks([offset + x for x in range(len(options.xticks))],
                     options.xticks,
                     rotation="vertical",
                     fontsize="8")
    else:
        if col_headers and len(col_headers) < 100:
            pylab.xticks([offset + x for x in range(len(col_headers))],
                         col_headers,
                         rotation="vertical",
                         fontsize="8")

    if options.yticks:
        pylab.yticks([offset + y for y in range(len(options.yticks))],
                     options.yticks,
                     fontsize="8")
    else:
        if row_headers and len(row_headers) < 100:
            pylab.yticks([offset + y for y in range(len(row_headers))],
                         row_headers,
                         fontsize="8")

 def savepng(pre, img, title=None, **kwargs):
     fn = '%s-%s.png' % (pre, idstr)
     print 'Saving', fn
     plt.clf()
     plt.imshow(img, **kwargs)
     ax = plt.axis()
     if debug:
         print len(xplotx),len(allobjx)
         for i,(objx,objy,objc) in enumerate(zip(allobjx,allobjy,allobjc)):
             plt.plot(objx,objy,'-',c=objc)
             tempx = []
             tempx.append(xplotx[i])
             tempx.append(objx[0])
             tempy = []
             tempy.append(xploty[i])
             tempy.append(objy[0])
             plt.plot(tempx,tempy,'-',c='purple')
         plt.plot(pointx,pointy,'y.')
         plt.plot(xplotx,xploty,'xg')
     plt.axis(ax)
     if title is not None:
         plt.title(title)
     plt.colorbar()
     plt.gray()
     plt.savefig(fn)

def psfplots():
	tpsf = wise.get_psf_model(1, pixpsf=True)
	
	psfp = tpsf.getPointSourcePatch(0, 0)
	psf = psfp.patch
	
	psf /= psf.sum()
	
	plt.clf()
	plt.imshow(np.log10(np.maximum(1e-5, psf)), interpolation='nearest', origin='lower')
	plt.colorbar()
	ps.savefig()
	
	h,w = psf.shape
	cx,cy = w/2, h/2
	
	X,Y = np.meshgrid(np.arange(w), np.arange(h))
	R = np.sqrt((X - cx)**2 + (Y - cy)**2)
	plt.clf()
	plt.semilogy(R.ravel(), psf.ravel(), 'b.')
	plt.xlabel('Radius (pixels)')
	plt.ylabel('PSF value')
	plt.ylim(1e-8, 1.)
	ps.savefig()
	
	plt.clf()
	plt.loglog(R.ravel(), psf.ravel(), 'b.')
	plt.xlabel('Radius (pixels)')
	plt.ylabel('PSF value')
	plt.ylim(1e-8, 1.)
	ps.savefig()
	
	print('PSF norm:', np.sqrt(np.sum(np.maximum(0, psf)**2)))
	print('PSF max:', psf.max())

def dispims_color(M, border=0, bordercolor=[0.0, 0.0, 0.0], savePath=None, *imshow_args, **imshow_keyargs):
    """ Display an array of rgb images. 

    The input array is assumed to have the shape numimages x numpixelsY x numpixelsX x 3
    """
    bordercolor = numpy.array(bordercolor)[None, None, :]
    numimages = len(M)
    M = M.copy()
    for i in range(M.shape[0]):
        M[i] -= M[i].flatten().min()
        M[i] /= M[i].flatten().max()
    height, width, three = M[0].shape
    assert three == 3
    
    n0 = numpy.int(numpy.ceil(numpy.sqrt(numimages)))
    n1 = numpy.int(numpy.ceil(numpy.sqrt(numimages)))
    im = numpy.array(bordercolor)*numpy.ones(
                             ((height+border)*n1+border,(width+border)*n0+border, 1),dtype='<f8')
    for i in range(n0):
        for j in range(n1):
            if i*n1+j < numimages:
                im[j*(height+border)+border:(j+1)*(height+border)+border,
                   i*(width+border)+border:(i+1)*(width+border)+border,:] = numpy.concatenate((
                  numpy.concatenate((M[i*n1+j,:,:,:],
                         bordercolor*numpy.ones((height,border,3),dtype=float)), 1),
                  bordercolor*numpy.ones((border,width+border,3),dtype=float)
                  ), 0)
    imshow_keyargs["interpolation"]="nearest"
    pylab.imshow(im, *imshow_args, **imshow_keyargs)
    
    if savePath == None:
        pylab.show()
    else:
        pylab.savefig(savePath)

 def img_create(self):
     '''
     Create the output jpg of the file.
     :return:
     '''
     pylab.imshow(self._dcm.pixel_array, cmap=pylab.cm.bone)
     pylab.savefig(self._str_outputImageFile)

    def _draw_solution(self):
        """ plot the current solution on our optional visualization """
        myg = self.grids[self.current_level].grid

        v = self.grids[self.current_level].get_var("v")

        pylab.imshow(
            numpy.transpose(v[myg.ilo : myg.ihi + 1, myg.jlo : myg.jhi + 1]),
            interpolation="nearest",
            origin="lower",
            extent=[self.xmin, self.xmax, self.ymin, self.ymax],
        )

        # pylab.xlabel("x")
        pylab.ylabel("y")

        if self.current_level == self.nlevels - 1:
            pylab.title(r"solving $L\phi = f$")
        else:
            pylab.title(r"solving $Le = r$")

        formatter = matplotlib.ticker.ScalarFormatter(useMathText=True)
        cb = pylab.colorbar(format=formatter, shrink=0.5)

        cb.ax.yaxis.offsetText.set_fontsize("small")
        cl = pylab.getp(cb.ax, "ymajorticklabels")
        pylab.setp(cl, fontsize="small")