def connection_field_plot_continuous(self,index,afferent=True,density=30):
     weights =  self.proj.getWeights(format='array')
     x = []
     y = []
     w = []
     
     if afferent:
            weights = weights[:,index].ravel()
            p = self.proj.pre
     else:
            weights = weights[index,:].ravel()
            p = self.proj.post
     
     for (ww,i) in zip(weights,numpy.arange(0,len(weights),1)):
               x.append(p.positions[0][i])
               y.append(p.positions[1][i])
               w.append(ww)
         
     bx = min([min(p.positions[0]),min(p.positions[0])])  
     by = max([max(p.positions[1]),max(p.positions[1])])  
     xi = numpy.linspace(min(p.positions[0]),max(p.positions[0]),100)
     yi = numpy.linspace(min(p.positions[1]),max(p.positions[1]),100)
     zi = griddata(x,y,w,xi,yi)
     pylab.figure()
     pylab.imshow(zi)
     pylab.title('Connection field from %s to %s of neuron %d' % (self.source.name,self.target.name,index))
     pylab.colorbar()

    def displayResults(self,res, cm=pylab.cm.gray, title='Specify a title'):
        if self.display:
		self.count=self.count+1
        	pylab.figure(self.count)
        	pylab.imshow(res, cm, interpolation='nearest')
        	pylab.colorbar()
        	pylab.title(title)

def getOptCandGamma(cv_train, cv_label):
    print "Finding optimal C and gamma for SVM with RBF Kernel"
    C_range = 10.0 ** np.arange(-2, 9)
    gamma_range = 10.0 ** np.arange(-5, 4)
    param_grid = dict(gamma=gamma_range, C=C_range)
    cv = StratifiedKFold(y=cv_label, n_folds=40)

    # Use the svm.SVC() as the cost function to evaluate parameter choices
    # NOTE: Perhaps we should run computations in parallel if needed. Does it
    # do that already within the class?
    grid = GridSearchCV(svm.SVC(), param_grid=param_grid, cv=cv)
    grid.fit(cv_train, cv_label)

    score_dict = grid.grid_scores_
    scores = [x[1] for x in score_dict]
    scores = np.array(scores).reshape(len(C_range), len(gamma_range))
    pl.figure(figsize=(8,6))
    pl.subplots_adjust(left=0.05, right=0.95, bottom=0.15, top=0.95)
    pl.imshow(scores, interpolation='nearest', cmap=pl.cm.spectral)
    pl.xlabel('gamma')
    pl.ylabel('C')
    pl.colorbar()
    pl.xticks(np.arange(len(gamma_range)), gamma_range, rotation=45)
    pl.yticks(np.arange(len(C_range)), C_range)
    pl.show()

    print "The best classifier is: ", grid.best_estimator_

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 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 plot_worker(jobq,mask,pid,lineshape,range):
    '''
    args[0] = array file name
    args[1] = output figure name
    if mask, where masked==0 is masked
    '''

    if lineshape:
        lines = shapefile.load_shape_list(lineshape)
    else:
        lines = None
    while True:
        #--get some args from the queue
        args = jobq.get()
        #--check if this is a sentenial
        if args == None:
            break
        #--load
        if args[2]:
            arr = np.fromfile(args[0],dtype=np.float32)
            arr.resize(bro.nrow,bro.ncol)
        else:
            arr = np.loadtxt(args[0])
        
        if mask != None:
            arr = np.ma.masked_where(mask==0,arr)        
        #print args[0],arr.min(),arr.max(),arr.mean()
        #--generic plotting
        fig = pylab.figure()
        ax = pylab.subplot(1,1,1,aspect='equal')
        
        if range:
            vmax = range[1]
            vmin = range[0]
        else:
            vmax = arr.max()
            vmin = arr.min()

        #p = ax.imshow(arr,interpolation='none')        
        p = ax.pcolor(bro.X,bro.Y,np.flipud(arr),vmax=vmax,vmin=vmin)
        pylab.colorbar(p)
        if lines:
            for line in lines:
                ax.plot(line[0,:],line[1,:],'k-',lw=1.0)
                #break
        ax.set_xticklabels([])
        ax.set_yticklabels([])
        ax.set_xlim(bro.plt_x)
        ax.set_ylim(bro.plt_y)
        ax.set_title(args[0])
        fmt = args[1].split('.')[-1]
        pylab.savefig(args[1],dpi=300,format=fmt)
        pylab.close(fig)
        #--mark this task as done
        jobq.task_done()
        print 'plot worker',pid,' finished',args[0]

    #--mark the sentenial as done
    jobq.task_done()
    return

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_x=%.4f,\  D_{-}=%.4f,\  D_{+}=%.4f,\ ch=%i$ ' %
                    (job.u_x, job.D_minus, job.D_plus, job.ch_objects))
        car = pylab.imshow(plane, vmin=-0.001, vmax=0.0, 
                           interpolation='nearest')
        pylab.contour(plane, contour_values, linestyles='dashed', 
                                             colors='white')
        pylab.grid(True)
        pylab.colorbar(car)
        #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 draw_heat_graph(getgraph, opts):
    # from pyevolve_graph script

    stage_points = getgraph()

    fg = pl.figure()
    ax = fg.add_subplot(111)

    pl.imshow(
        stage_points, aspect="auto", interpolation="gaussian",
        cmap=matplotlib.cm.__dict__["jet"])
    pl.title("Population scores along the generations")

    def labelfmt(x, pos=0):
        # there is surely a better way to do that
        return (float(x) == int(x)) and '%d' % (x) or ''
    ax.xaxis.set_major_formatter(pl.FuncFormatter(labelfmt))

    pl.xlabel('Generations -->')
    pl.ylabel('Sorted Population Results')
    pl.grid(True)
    pl.colorbar()

    if opts.outfile:
        fg.savefig(opts.outfile)

    if opts.show:
        pl.show()

def plot_C_gamma_grid_search(grid, C_range, gamma_range, score):
    '''
    Plots the scores computed on a grid. 
    
    Arguments: 
        grid - the grid search object created using GridSearchCV()
        C_range - the C parameter range 
        gamma_range - the gamma parameter range 
        score - the scoring function  
        
    
    '''

    # grid_scores_ contains parameter settings and scores
    # We extract just the scores
    scores = [x[1] for x in grid.grid_scores_]
    scores = np.array(scores).reshape(len(C_range), len(gamma_range))
    
    # draw heatmap of accuracy as a function of gamma and C
    pl.figure(figsize=(8, 6))
    pl.subplots_adjust(left=0.05, right=0.95, bottom=0.15, top=0.95)
    pl.imshow(scores, interpolation='nearest', cmap=pl.cm.spectral)
    pl.title("Grid search on C and gamma for best %s" % score)
    pl.xlabel('gamma')
    pl.ylabel('C')
    pl.colorbar()
    pl.xticks(np.arange(len(gamma_range)), gamma_range, rotation=45)
    pl.yticks(np.arange(len(C_range)), C_range)
    
    pl.show()

 def plothistory(self):
     a=self.a
     b=self.b
     plt.figure(figsize=(12,6))
     I=np.concatenate([a.T,np.array(np.nansum(a[:,:3],1),ndmin=2),
         np.array(np.nansum(a[:,3:6],1),ndmin=2),np.array(np.nansum(a[:,6:],1),ndmin=2)],axis=0)
     plt.plot(range(b.size),b,'rx',ms=8,mew=2)
     plt.plot([10.5,10.5],[-1,I.shape[1]],'r',lw=2)
     plt.imshow(I,interpolation='nearest',cmap='winter')
     plt.colorbar()
     ax=plt.gca()
     ax.set_yticks(range(I.shape[0]))
     ax.set_yticklabels(['']*a.shape[0]+['color','rel len','abs len'])
     c1=plt.Circle((-1.5,0),radius=0.4,color='blue',clip_on=False)
     c2=plt.Circle((-1.5,1),radius=0.4,color='white',clip_on=False)
     c3=plt.Circle((-1.5,2),radius=0.4,color='yellow',clip_on=False)
     ax.add_patch(c1);ax.add_patch(c2);ax.add_patch(c3);
     c1=plt.Rectangle((-2,3),1,0.2,color='white',clip_on=False)
     c2=plt.Rectangle((-2.5,4),1.5,0.2,color='white',clip_on=False)
     c3=plt.Rectangle((-3,5),2,0.2,color='white',clip_on=False)
     ax.add_patch(c1);ax.add_patch(c2);ax.add_patch(c3);
     c1=plt.Rectangle((-2,6),1,0.2,color='gray',clip_on=False)
     c2=plt.Rectangle((-2.5,7),1.5,0.2,color='gray',clip_on=False)
     c3=plt.Rectangle((-3,8),2,0.2,color='gray',clip_on=False)
     c4=plt.Rectangle((-3.5,9),2.5,0.2,color='gray',clip_on=False)
     ax.add_patch(c1);ax.add_patch(c2);ax.add_patch(c3);ax.add_patch(c4);
     print I[-3,-1]

def plot_samples_distance(dataset, sortbyattr=None):
    """Plot the euclidean distances between all samples of a dataset.

    Parameters
    ----------
    dataset : Dataset
      Providing the samples.
    sortbyattr : None or str
      If None, the samples distances will be in the same order as their
      appearance in the dataset. Alternatively, the name of a samples
      attribute can be given, which wil then be used to sort/group the
      samples, e.g. to investigate the similarity samples by label or by
      chunks.
    """
    if sortbyattr is not None:
        slicer = []
        for attr in dataset.sa[sortbyattr].unique:
            slicer += \
                get_samples_by_attr(dataset, sortbyattr, attr).tolist()
        samples = dataset.samples[slicer]
    else:
        samples = dataset.samples

    ed = np.sqrt(squared_euclidean_distance(samples))

    pl.imshow(ed, interpolation='nearest')
    pl.colorbar()

def plot_mtx(mtx=None, title=None, newfig=False, cbar=True, **kwargs):
    """
    ::

        static method for plotting a matrix as a time-frequency distribution (audio features)
    """
    if mtx is None or type(mtx) != np.ndarray:
        raise ValueError('First argument, mtx, must be a array')
    if newfig: P.figure()
    dbscale = kwargs.pop('dbscale', False) 
    bels = kwargs.pop('bels',False)
    norm = kwargs.pop('norm',False)
    normalize = kwargs.pop('normalize',False)
    origin=kwargs.pop('origin','lower')
    aspect=kwargs.pop('aspect','auto')
    interpolation=kwargs.pop('interpolation','nearest')
    cmap=kwargs.pop('cmap',P.cm.gray_r)
    clip=-100.
    X = scale_mtx(mtx, normalize=normalize, dbscale=dbscale, norm=norm, bels=bels)
    i_min, i_max = np.where(X.mean(1))[0][[0,-1]]
    X = X[i_min:i_max+1].copy()
    if dbscale or bels:
        if bels: clip/=10.
        P.imshow(P.clip(X,clip,0),origin=origin, aspect=aspect, interpolation=interpolation, cmap=cmap, **kwargs)
    else:
        P.imshow(X,origin=origin, aspect=aspect, interpolation=interpolation, cmap=cmap, **kwargs)
    if title:
        P.title(title,fontsize=16)
    if cbar:
        P.colorbar()
    P.yticks(np.arange(0,i_max+1-i_min,3),pc_labels[i_min:i_max+1:3],fontsize=14)
    P.xlabel('Tactus', fontsize=14)
    P.ylabel('MIDI Pitch', fontsize=14)
    P.grid()

def plot_dendrogram_and_matrix(linkage, matrix, color_threshold=None):
    # Compute and plot dendrogram.
    fig = pylab.figure(figsize=(20,20))
    axdendro = fig.add_axes([0.09,0.1,0.2,0.8])
    dendrogram = sch.dendrogram(linkage, color_threshold=color_threshold, orientation='right')
    axdendro.set_xticks([])
    axdendro.set_yticks([])

    # Plot distance matrix.
    axmatrix = fig.add_axes([0.3,0.1,0.6,0.8])
    index = dendrogram['leaves']
    D = matrix[:]
    D = D[index,:]
    D = D[:,index]
    im = axmatrix.matshow(D, aspect='auto', origin='lower')
    axmatrix.set_xticks([])
    axmatrix.set_yticks([])

    # Plot colorbar.
    axcolor = fig.add_axes([0.91,0.1,0.02,0.8])
    pylab.colorbar(im, cax=axcolor)

    # Display and save figure.
    fig.show()
    raw_input()

    return dendrogram

def plot_vel_vs_h3(maps, out_suffix=''):
    xaxis = (np.arange(len(maps.velocity[0])) - ppxf_m31.bhpos_pix[0]) * 0.05
    yaxis = (np.arange(len(maps.velocity)) - ppxf_m31.bhpos_pix[1]) * 0.05
    yy, xx = np.meshgrid(xaxis, yaxis)
    radius = np.hypot(xx, yy)
    good = np.where((np.abs(yy) < 0.5) & (np.abs(xx) < 1.0))
        
    plt.scatter(maps.velocity[good], maps.h3[good], c=maps.sigma[good], s=5,
                    marker='o', vmin=0, vmax=450)
    plt.xlim(-700, 0)
    plt.ylim(-0.5, 0.5)
    plt.colorbar(label='Sigma (km/s)')
    plt.axhline(linestyle='--', color='grey')
    plt.xlabel('Velocity (km/s)')
    plt.ylabel('h3')
    plt.savefig(plot_dir + 'vel_vs_h3' + out_suffix + '.png')

    plt.clf()
    plt.scatter(maps.sigma[good], maps.h3[good], c=maps.velocity[good], s=5,
                    marker='o', vmin=-700, vmax=0)
    plt.xlim(0, 450)
    plt.ylim(-0.5, 0.5)
    plt.colorbar(label='Velocity (km/s)')
    plt.axhline(linestyle='--', color='grey')
    plt.xlabel('Sigma (km/s)')
    plt.ylabel('h3')
    plt.savefig(plot_dir + 'sig_vs_h3' + out_suffix + '.png')

    return

    def plot_plasma(self):

        P.tricontourf(self.rzt[:, 0], self.rzt[:, 1],
                      self.tris, self.beta, 1001, zorder=0)
        cticks = P.linspace(0.0, 0.2, 5)
        P.colorbar(ticks=cticks, format='%.2f')
        P.jet()

def plot_box_data(field,redshift):
    nf1 = get_box_data(field,redshift)
    chosenIndex=200
    nf1 = nf1.reshape((400,400,400))[chosenIndex,:,:]
    plt.imshow(nf1)
    plt.colorbar()
    plt.show()

def pressx(ifile, varkey, options, before = '', after = ''):
    import pylab as pl
    from matplotlib.colors import Normalize, LogNorm
    outpath = getattr(options, 'outpath', '.')
    vert = getpresbnds(ifile)
    var = ifile.variables[varkey]
    dims = [(k, l) for l, k in zip(var[:].shape, var.dimensions) if l > 1]
    if len(dims) > 2:
        raise ValueError('Press-x can have 2 non-unity dimensions; got %d - %s' % (len(dims), str(dims)))
    if options.logscale:
        norm = LogNorm()
    else:
        norm = Normalize()
    exec(before)
    ax = pl.gca()
    print(varkey, end = '')
    vals = var[:].squeeze()
    x = np.arange(vals.shape[1])
    patches = ax.pcolor(x, vert, vals, norm = norm)
    #ax.set_xlabel(X.units.strip())
    #ax.set_ylabel(Y.units.strip())
    pl.colorbar(patches)
    ax.set_ylim(vert.max(), vert.min())
    ax.set_xlim(x.min(), x.max())
    fmt = 'png'
    figpath = os.path.join(outpath + '_PRESX_' + varkey + '.' + fmt)
    exec(after)
    pl.savefig(figpath)
    print('Saved fig', figpath)
    return figpath

def main():
    base_path = "/caps2/tsupinie/1kmf-control/"
    temp = goshen_1km_temporal(start=14400, end=14400)
    grid = goshen_1km_grid()
    n_ens_members = 40

    np.seterr(all='ignore')

    ens = loadEnsemble(base_path, [ 11 ], temp.getTimes(), ([ 'pt', 'p' ], computeDensity))
    ens = ens[0, 0]

    zs = decompressVariable(nio.open_file("%s/ena001.hdfgrdbas" % base_path, mode='r', format='hdf').variables['zp'])
    xs, ys = grid.getXY()
    xs = xs[np.newaxis, ...].repeat(zs.shape[0], axis=0)
    ys = ys[np.newaxis, ...].repeat(zs.shape[0], axis=0)

    eff_buoy = effectiveBuoyancy(ens, (zs, ys, xs), plane={'z':10})
    print eff_buoy

    pylab.figure()
    pylab.contourf(xs[0], ys[0], eff_buoy[0], cmap=matplotlib.cm.get_cmap('RdBu_r'))
    pylab.colorbar()

    grid.drawPolitical()

    pylab.suptitle("Effective Buoyancy")
    pylab.savefig("eff_buoy.png")
    pylab.close()
    return

 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 VisualizeAlm(alm,figno=1,max_l=None):
    """ Visualize a healpy a_lm vector """
    lmax = hp.Alm.getlmax(f_lm.size)
    l,m = hp.Alm.getlm(lmax)
    mag = np.zeros([lmax+1,lmax+1])
    phs = np.zeros([lmax+1,lmax+1])
    mag[m,l] = np.abs(alm)
    phs[m,l] = np.angle(alm)
    cl = hp.alm2cl(alm)
    # Decide the range of l to plot
    if max_l != None:
        max_l = (max_l if (max_l <= lmax) else lmax)
    else:
        max_l = lmax 
    print max_l
    plt.figure(figno)
    plt.clf()
    plt.subplot(211)
    plt.imshow(mag[0:max_l,0:max_l],interpolation='nearest',origin='lower')
    plt.colorbar()
    plt.subplot(212)
    plt.imshow(phs[0:max_l,0:max_l],interpolation='nearest',origin='lower')
    plt.colorbar()
    # plt.subplot(313)
    #plt.semilogy(cl[0:max_l])
    return {'mag':mag,'phs':phs,'cl':cl}