def plotGeometry(self,SNP,PCAD,title):
        fig=plt.figure(figsize=self.figsize, dpi=self.dpi);plt.ioff()
        SNP.loc['Reference0']=np.zeros(SNP.shape[1],dtype=int)
        plt.subplot(1,2,1)
        D=pd.DataFrame(None,index=SNP.index,columns=SNP.index)
        for i in SNP.index:
            for j in SNP.index:
                D.loc[i,j]= sum(np.logical_xor(SNP.loc[i],SNP.loc[j]))
        D=D.astype(float)
        im=plt.imshow(D,interpolation='nearest',cmap='Reds')
        plt.gca().xaxis.tick_top()
        x=np.arange(D.index.shape[0])
        plt.yticks(x,map(lambda x: x.replace('mdio','') ,D.index.values))
        plt.xticks(x,map(lambda x: x.replace('mdio','') ,D.columns))
        plt.colorbar(im)
        plt.gca().tick_params(axis='both', which='major', labelsize=10)
        plt.title('Pairwise Hamming Distance',y=1.03)

        
        plt.subplot(1,2,0)
        D=PCAD.astype(float)
        im=plt.imshow(D,interpolation='nearest',cmap='Reds')
        plt.gca().xaxis.tick_top()
        x=np.arange(D.index.shape[0])
        plt.yticks(x,map(lambda x: x.replace('mdio','') ,D.index.values))
        plt.xticks(x,map(lambda x: x.replace('mdio','') ,D.columns))
        plt.colorbar(im)
        plt.gca().tick_params(axis='both', which='major', labelsize=10)
        plt.title('Euclidean Distance in PC3',y=1.03)
        plt.suptitle('Figure {}. {}'.format(self.fignumber, title),fontsize=self.titleSizeSup); self.pdf.savefig(fig);self.fignumber+=1

    def rmsdSpreadSubplot(multiplier=1.0, layout=(-1, -1)):
        rmsd_data   = dict( (e, rad_data[e]['innov'][quant])  for e in rad_data.iterkeys() )
        spread_data = dict( (e, rad_data[e]['spread'][quant]) for e in rad_data.iterkeys() )

        times = temp.getTimes()
        n_t = len(times)

        for exp, exp_name in exp_names.iteritems():
            pylab.plot(sawtooth(times, times)[:(n_t + 1)], rmsd_data[exp][:(n_t + 1)], color=colors[exp], linestyle='-')
            pylab.plot(times[(n_t / 2):], rmsd_data[exp][n_t::2], color=colors[exp], linestyle='-')
 
        for exp, exp_name in exp_names.iteritems():
            pylab.plot(sawtooth(times, times)[:(n_t + 1)], spread_data[exp][:(n_t + 1)], color=colors[exp], linestyle='--')
            pylab.plot(times[(n_t / 2):], spread_data[exp][n_t::2], color=colors[exp], linestyle='--')

        ylim = pylab.ylim()
        pylab.plot(times, -1 * np.ones((len(times),)), color='#999999', linestyle='-', label="RMS Innovation")
        pylab.plot(times, -1 * np.ones((len(times),)), color='#999999', linestyle='--', label="Spread")

        pylab.axhline(y=7, color='k', linestyle=':')
        pylab.axvline(x=14400, color='k', linestyle=':')

        pylab.ylabel("RMS Innovation/Spread (dBZ)", size='large')

        pylab.xlim(times[0], times[-1])
        pylab.ylim(ylim)

        pylab.legend(loc=4)

        pylab.xticks(times[::2], [ "" for t in times[::2] ])
        pylab.yticks(size='x-large')
        return

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 RosenbrockTest():
    nDim = 3
    numOfParticles = 20
    maxIteration = 200
    minX = array([-5.0]*nDim)
    maxX = array([5.0]*nDim)
    maxV = 0.2*(maxX - minX)
    minV = -1.0*maxV
    numOfTrial = 20
    gBest = array([0.0]*maxIteration)
    for i in xrange(numOfTrial):
        p1 = RPSO.PSOProblem(nDim, numOfParticles, maxIteration, minX, maxX, minV, maxV, RPSO.Rosenbrock)
        p1.run()
        gBest = gBest + p1.gBestArray[:maxIteration]
    gBest = gBest / numOfTrial
    pylab.title('$G_{best}$ over 20 trials')
    pylab.xlabel('The $N^{th}$ Iteratioin')
    pylab.ylabel('Average gBest over '+str(numOfTrial)+' runs (logscale)')
    pylab.grid(True)
#    pylab.yscale('log')
    ymin, ymax = -1.5, 2.5
    ystep = 0.5
    pylab.ylim(ymin, ymax)
    yticks = linspace(ymin, ymax, (ymax-ymin)/ystep+1)
    pylab.yticks(tuple(yticks),tuple(map(str,yticks)))
    pylab.plot(range(maxIteration), log10(gBest),'-', label='Global best')
    pylab.legend()
    pylab.show()

def command(args):
    from pylab import bar, yticks, subplots_adjust, show
    from numpy import arange

    import sr.tools.bom.bom as bom
    import sr.tools.bom.parts_db as parts_db

    db = parts_db.get_db()
    m = bom.MultiBoardBom(db)
    m.load_boards_args(args.arg)
    m.prime_cache()

    prices = []

    for srcode, pg in m.items():
        if srcode == "sr-nothing":
            continue

        prices.append((srcode, pg.get_price()))

    prices.sort(key=lambda x: x[1])

    bar(0, 0.8, bottom=range(0, len(prices)), width=[x[1] for x in prices],
        orientation='horizontal')

    yticks(arange(0, len(prices)) + 0.4, [x[0] for x in prices])

    subplots_adjust(left=0.35)

    show()

def tracks_movie(base, skip=1, frames=500, size=10):
    """
    A movie of each particle as a point
    """
    conf, track, pegs = load(base)

    fig = pl.figure(figsize=(size,size*conf['top']/conf['wall']))
    plot = None

    for t in xrange(1,max(frames, track.shape[1]/skip)):
        tmp = track[:,t*skip,:]
        if not ((tmp[:,0] > 0) & (tmp[:,1] > 0) & (tmp[:,0] < conf['wall']) & (tmp[:,1] < conf['top'])).any():
            continue

        if plot is None:
            plot = pl.plot(tmp[:,0], tmp[:,1], 'k,', alpha=1.0, ms=0.1)[0]
            pl.xticks([])
            pl.yticks([])
            pl.xlim(0,conf['wall'])
            pl.ylim(0,conf['top'])
            pl.tight_layout()
        else:
            plot.set_xdata(tmp[:,0])
            plot.set_ydata(tmp[:,1])
        pl.draw()
        pl.savefig(base+'-movie-%05d.png' % (t-1))

def plot_importances(imp, clfName, obj):
    imp=np.vstack(imp)
    print imp
    mean_importance = np.mean(imp,axis=0)
    std_importance = np.std(imp,axis=0)
    indices = np.argsort(mean_importance)[::-1]
    print indices
    print featureNames
    featureList = []
    # num_features = len(featureNames)
    print("Feature ranking:")
    for f in range(num_features):
        featureList.append(featureNames[indices[f]])
        print("%d. feature %s (%.2f)" % (f, featureNames[indices[f]], mean_importance[indices[f]]))
    fig = pl.figure(figsize=(8,6),dpi=150)
    pl.title("Feature importances",fontsize=30)
    pl.bar(range(num_features), mean_importance[indices],
            yerr = std_importance[indices], color=paired[0], align="center",
            edgecolor=paired[0],ecolor=paired[1])
    pl.xticks(range(num_features), featureList, size=15,rotation=90)
    pl.ylabel("Importance",size=30)
    pl.yticks(size=20)
    pl.xlim([-1, num_features])
    # fix_axes()
    pl.tight_layout()
    save_path = 'plots/'+obj+'/'+clfName+'_feature_importances.pdf'
    fig.savefig(save_path)

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_multiple_roc(rocList,title='',labels=None, include_baseline=False, equal_aspect=True):
	""" Plots multiple ROC curves on the same chart. 
		Parameters:
			rocList: the list of ROCData objects
			title: The tile of the chart
			labels: The labels of each ROC curve
			include_baseline: if it's  True include the random baseline
			equal_aspect: keep equal aspect for all roc curves
	"""
	pylab.clf()
	pylab.ylim((0,1))
	pylab.xlim((0,1))
	pylab.xticks(pylab.arange(0,1.1,.1))
	pylab.yticks(pylab.arange(0,1.1,.1))
	pylab.grid(True)
	if equal_aspect:
		cax = pylab.gca()
		cax.set_aspect('equal')
	pylab.xlabel("1 - Specificity")
	pylab.ylabel("Sensitivity")
	pylab.title(title)
	if not labels:
		labels = [ '' for x in rocList]
	_remove_duplicate_styles(rocList)
	for ix, r in enumerate(rocList):
		pylab.plot([x[0] for x in r.derived_points], [y[1] for y in r.derived_points], r.linestyle, linewidth=1, label=labels[ix])
	if include_baseline:
		pylab.plot([0.0,1.0], [0.0, 1.0], 'k-', label= 'random')
	if labels:
		pylab.legend(loc='lower right')
		
	pylab.show()

    def gui_repr(self):
        """Generate a GUI to represent the sentence alignments
        """
        if __pylab_loaded__:
            fig_width = max(len(self.text_e), len(self.text_f)) + 1
            fig_height = 3
            pylab.figure(figsize=(fig_width*0.8, fig_height*0.8), facecolor='w')
            pylab.box(on=False)
            pylab.subplots_adjust(left=0, right=1, bottom=0, top=1)
            pylab.xlim(-1, fig_width - 1)
            pylab.ylim(0, fig_height)
            pylab.xticks([])
            pylab.yticks([])

            e = [0 for _ in xrange(len(self.text_e))]
            f = [0 for _ in xrange(len(self.text_f))]
            for (i, j) in self.align:
                e[i] = 1
                f[j] = 1
                # draw the middle line
                pylab.arrow(i, 2, j - i, -1, color='r')
            for i in xrange(len(e)):
                # draw e side line
                pylab.text(i, 2.5, self.text_e[i], ha = 'center', va = 'center',
                        rotation=30)
                if e[i] == 1:
                    pylab.arrow(i, 2.5, 0, -0.5, color='r', alpha=0.3, lw=2)
            for i in xrange(len(f)):
                # draw f side line
                 pylab.text(i, 0.5, self.text_f[i], ha = 'center', va = 'center',
                        rotation=30)
                 if f[i] == 1:
                    pylab.arrow(i, 0.5, 0, 0.5, color='r', alpha=0.3, lw=2)

            pylab.draw()

def animation_compare(filename1, filename2, prefix, tend):
    num = 0
    for t in timestep(0,tend):
        t1,s1 = FieldInitializer.LoadState(filename1, t)
        t2,s2 = FieldInitializer.LoadState(filename2, t)
        rot1 = s1.CalculateRotationRodrigues()['z']
        rot2 = s2.CalculateRotationRodrigues()['z']

        pylab.figure(figsize=(5.12*2,6.12))
        rho = s2.CalculateRhoFourier().modulus()
        rot = s2.CalculateRotationRodrigues()['z']
        pylab.subplot(121)
        #pylab.imshow(rho*(rho>0.5), alpha=1, cmap=pylab.cm.bone_r)
        pylab.imshow(rot1)
        pylab.xticks([])
        pylab.yticks([])
        pylab.xlabel(r'$d_0$', fontsize=25)
        pylab.subplot(122)
        pylab.imshow(rot2)
        pylab.xticks([])
        pylab.yticks([])
        pylab.subplots_adjust(0.,0.,1.,1.,0.01,0.05)
        pylab.xlabel(r'$d_0/2$', fontsize=25)
        pylab.suptitle("Nabarro-Herring", fontsize=25)
        pylab.savefig("%s%04i.png" %(prefix, num))
        pylab.close('all')
        num = num + 1

    def default_frame21(self, data_label, x_tpl, y1_tpl,
                        y2_tpl, pltstyle_dict):
        x, labelx = x_tpl
        y1, label1 = y1_tpl
        y2, label2 = y2_tpl

        self.ax1 = pylab.subplot(211)
        self.ax1.plot(x, y1, label=data_label, **pltstyle_dict)
        pylab.ylabel(label1, fontsize=20, horizontalalignment='left')
        pylab.yticks(fontsize=13)
        self.ax1.yaxis.set_label_coords(-.12, 0.5)
        self.ax1.legend()

        self.ax2 = pylab.subplot(212)
        self.ax2.plot(x, y2, **pltstyle_dict)
        #self.ax2.yaxis.tick_right()
        #pylab.ylim(0,3)
        #self.ax2.yaxis.set_label_coords(1.12, 0.5)
        pylab.ylabel(label2, fontsize=20)

        #xticklabels = self.ax1.get_xticklabels()+self.ax2.get_xticklabels()
        #pylab.setp(xticklabels, visible=False)
        #pylab.ylabel(label3, fontsize = 20)
        pylab.xlabel(labelx, fontsize=20)
        pylab.xticks(fontsize=13)

    def default_frame31shareX(self, data_label, x_tpl, y1_tpl,
                              y2_tpl, y3_tpl, pltstyle_dict):
        x, labelx = x_tpl
        y1, label1 = y1_tpl
        y2, label2 = y2_tpl
        y3, label3 = y3_tpl

        pylab.subplots_adjust(hspace=0.001)
        self.ax1 = pylab.subplot(311)
        self.ax1.plot(x, y1, label=data_label, **pltstyle_dict)
        pylab.ylabel(label1, fontsize=20, horizontalalignment='left')
        pylab.yticks(fontsize=13)
        #self.ax1.yaxis.set_label_coords(-.12, 0.5)
        self.ax1.yaxis.set_label_coords(-.10, 0.25)

        pylab.ylim(-30, -10)
        self.ax1.legend()

        self.ax2 = pylab.subplot(312, sharex=self.ax1)
        self.ax2.plot(x, y2, **pltstyle_dict)
        self.ax2.yaxis.tick_right()
        pylab.ylim(0, 3)
        self.ax2.yaxis.set_label_coords(1.12, 0.5)
        pylab.ylabel(label2, fontsize=20)

        self.ax3 = pylab.subplot(313, sharex=self.ax1)
        self.ax3.plot(x, y3, **pltstyle_dict)

        xticklabels = self.ax1.get_xticklabels()+self.ax2.get_xticklabels()
        pylab.setp(xticklabels, visible=False)
        pylab.ylabel(label3, fontsize=20)
        pylab.xlabel(labelx, fontsize=20)
        pylab.xticks(fontsize=13)

def decorate(mean):
    pl.legend(loc='upper center', bbox_to_anchor=(.5,-.3))
    xmin, xmax, ymin, ymax = pl.axis()
    pl.vlines([mean], -ymax, ymax*10, linestyle='dashed', zorder=20)
    pl.xticks([0, .5, 1])
    pl.yticks([])
    pl.axis([-.01, 1.01, -ymax*.05, ymax*1.01])

def plot_question(fname, question_text, data):
    import pylab
    import numpy as np
    from matplotlib.font_manager import FontProperties
    from matplotlib.text import Text
    pylab.figure().clear()
    pylab.title(question_text)
    #pylab.xlabel("Verteilung")
    #pylab.subplot(101)
    if True or len(data) < 3:
        width = 0.95
        pylab.bar(range(len(data)), [max(y, 0.01) for x, y in data], 0.95, color="g")
        pylab.xticks([i+0.5*width for i in range(len(data))], [x for x, y in data])
        pylab.yticks([0, 10, 20, 30, 40, 50])
        #ind = np.arange(len(data))
        #pylab.bar(ind, [y for x, y in data], 0.95, color="g")
        #pylab.ylabel("#")
        #pylab.ylim(ymax=45)
        #pylab.ylabel("Antworten")
        #pylab.xticks(ind+0.5, histo.get_ticks())
        #pylab.legend(loc=3, prop=FontProperties(size="smaller"))
        ##pylab.grid(True)
    else:
        pylab.pie([max(y, 0.1) for x, y in data], labels=[x for x, y in data], autopct="%.0f%%")
    pylab.savefig(fname, format="png", dpi=75)

	def plot(self,title='',include_baseline=False,equal_aspect=True):
		""" Method that generates a plot of the ROC curve 
			Parameters:
				title: Title of the chart
				include_baseline: Add the baseline plot line if it's True
				equal_aspect: Aspects to be equal for all plot
		"""
		
		pylab.clf()
		pylab.plot([x[0] for x in self.derived_points], [y[1] for y in self.derived_points], self.linestyle)
		if include_baseline:
			pylab.plot([0.0,1.0], [0.0,1.0],'k-.')
		pylab.ylim((0,1))
		pylab.xlim((0,1))
		pylab.xticks(pylab.arange(0,1.1,.1))
		pylab.yticks(pylab.arange(0,1.1,.1))
		pylab.grid(True)
		if equal_aspect:
			cax = pylab.gca()
			cax.set_aspect('equal')
		pylab.xlabel('1 - Specificity')
		pylab.ylabel('Sensitivity')
		pylab.title(title)
		
		pylab.show()

def plot_file_color(base, thin=True, start=0, size=14, save=False):
    conf, track, pegs = load(base)

    fig = pl.figure(figsize=(size,size*conf['top']/conf['wall']))

    track = track[start:]
    x = track[:,0];   y = track[:,1]
    t = np.linspace(0,1,x.shape[0])
    points = np.array([x,y]).transpose().reshape(-1,1,2)
    segs = np.concatenate([points[:-1],points[1:]],axis=1)
    lc = LineCollection(segs, linewidths=0.25, cmap=pl.cm.coolwarm)
    lc.set_array(t)
    pl.gca().add_collection(lc)

    #pl.scatter(x, y, c=np.arange(len(x)),linestyle='-',cmap=pl.cm.coolwarm)
    #pl.plot(track[-1000000:,0], track[-1000000:,1], '-', linewidth=0.0125, alpha=0.8)
    for peg in pegs:
        pl.gca().add_artist(pl.Circle(peg, conf['radius'], color='k', alpha=0.3))
    pl.xlim(0, conf['wall'])
    pl.ylim(0, conf['top'])
    pl.xticks([])
    pl.yticks([])
    pl.tight_layout()
    pl.show()
    if save:
        pl.savefig(base+".png", dpi=200)

def plotLDDecaySelection3d(ax, sweep=False):
    import pylab as plt; import matplotlib as mpl;mpl.rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size':16}) ;    mpl.rc('text', usetex=True)

    def neutral(ld0, t, d, r=2 * 1e-8):
        if abs(d) <= 5e3:
            d = np.sign(d) * 5e3
        if d == 0:
            d = 5e3
        return ((np.exp(-2 * r * t * abs(d)))) * ld0

    t = np.arange(0, 200 + 1., 2)
    L=1e6+1
    pos=500000
    r=2*1e-8
    ld0 = 0.5
    s = 0.05
    nu0 = 0.1
    positions=np.arange(0,L,1000)
    dist=(positions - pos)
    T, D = np.meshgrid(t, dist)
    if not sweep:
        zs = np.array([neutral(ld0, t, d) for t, d in zip(np.ravel(T), np.ravel(D))])
    else:
        zs = np.array([LD(t, ld0, s, nu0, r, abs(d), 0) for t, d in zip(np.ravel(T), np.ravel(D))])
    Z = zs.reshape(T.shape)
    ax.plot_surface(T, D, Z,cmap=mpl.cm.autumn)
    ax.set_xlabel('Generations')
    ax.set_ylabel('Position')
    plt.yticks(plt.yticks()[0][1:-1],map(lambda x:'{:.0f}K'.format((pos+(x))/1000),plt.yticks()[0][1:-1]))
    plt.ylim([-500000,500000])
    ax.set_zlabel(r"$|\rho_t|$")
    pplt.setSize(plt.gca(), fontsize=6)
    plt.axis('tight');

def task1():
    '''
    Task 1
    Generate, transform and plot gaussian data
    '''
    X = generate_data(100)
    X2 = scale_data(X)
    X3 = standardise_data(X)

    # Plot data
    # Your code here
    # Hint: Use the functions pl.scatter(x[0,:],[1,:],c='r'), pl.hold(True),
    # pl.legend, pl.title, pl.xlabel, pl.ylabel
    fig = pl.figure()
    ax1 = fig.add_subplot(111)
    print ax1.scatter(X[0], X[1], c='y', label='Raw data')
    ax1.scatter(X2[0], X2[1], c='r', label='Scaled data')
    ax1.scatter(X3[0], X3[1], c='b', label='Standardised data')
    pl.title('Simple transformations of Gaussian Data')
    pl.xlabel('x')
    pl.ylabel('y')
    pl.xticks(range(-6,10,2))
    pl.yticks(range(-4,5,1))
    # ax1.title('Simple transformations of Gaussian Data') ???? -> does not work
    ax1.legend(scatterpoints=1)
    pl.savefig('task_1_340528_341455.pdf')

    def plot_bases(self, autoscale=True, stampsize=None):
        import pylab as plt

        N = len(self.psfbases)
        cols = int(np.ceil(np.sqrt(N)))
        rows = int(np.ceil(N / float(cols)))
        plt.clf()
        plt.subplots_adjust(hspace=0, wspace=0)

        cut = 0
        if stampsize is not None:
            H, W = self.shape
            assert H == W
            cut = max(0, (H - stampsize) / 2)

        ima = dict(interpolation="nearest", origin="lower")
        if autoscale:
            mx = self.psfbases.max()
            ima.update(vmin=-mx, vmax=mx)
        nil, xpows, ypows = self.polynomials(0.0, 0.0, powers=True)
        for i, (xp, yp, b) in enumerate(zip(xpows, ypows, self.psfbases)):
            plt.subplot(rows, cols, i + 1)

            if cut > 0:
                b = b[cut:-cut, cut:-cut]
            if autoscale:
                plt.imshow(b, **ima)
            else:
                mx = np.abs(b).max()
                plt.imshow(b, vmin=-mx, vmax=mx, **ima)
            plt.xticks([])
            plt.yticks([])
            plt.title("x^%i y^%i" % (xp, yp))
        plt.suptitle("PsfEx eigen-bases")