def regressionANN(mypathforResults,predicted,observed,regress,variable_to_fill, Site_ID,units,list_out,index_str):    
    for index, item in enumerate(list_out):
	fig=pl.figure(4, figsize=(16, 12), dpi=80, facecolor='w', edgecolor='k')
	ANN_label=str(item+"_NN") 
        graphtext1=str('slope      ' + str("{0:.2f}".format(regress[index][0])) +'\n' +
                       'intercept  ' + str("{0:.2f}".format(regress[index][1])) +'\n' +
                       'r-value    ' + str("{0:.2f}".format(regress[index][2])) +'\n' +
                       'p-value    ' + str("{0:.2f}".format(regress[index][3])) +'\n' +
                       'slope SE   ' + str("{0:.2f}".format(regress[index][4])) +'\n' +
                       'estim. SE  ' + str("{0:.2f}".format(regress[index][5]))         )  
        pl.figtext(0.7,0.6,graphtext1, bbox=dict())
        pl.plot(observed[:,index], predicted[:,index], 'o', label='targets vs. outputs')    
        slope = regress[index][0]; intercept = regress[index][1]
    
        x = np.linspace(min(observed[:,index]),max(observed[:,index]))
        y = slope * x + intercept
        pl.plot(x, y, linewidth = 2, label = 'regression line')
        pl.legend()
        pl.title('Tower vs ANN for '+item+' at ' +Site_ID+ ' index '+index_str)
        pl.xlabel('Tower ' + '('+units+')')
        pl.ylabel('ANN ' + '('+units+')')
        pl.legend(shadow=True, fancybox=True,loc='best')
        
        pl.savefig(mypathforResults+'/'+'Tower vs ANN for '+item+' at ' +Site_ID+ ' index '+index_str)
        #pl.show()
        pl.close(4)    
	time.sleep(2)

    def OnCalcShift(self, event):
        if (len(self.PSFLocs) > 0):
            import pylab
            
            x,y,z = self.PSFLocs[0]
            
            z_ = numpy.arange(self.image.data.shape[2])*self.image.mdh['voxelsize.z']*1.e3
            z_ -= z_.mean()
            
            pylab.figure()
            p_0 = 1.0*self.image.data[x,y,:,0].squeeze()
            p_0  -= p_0.min()
            p_0 /= p_0.max()

            #print (p_0*z_).sum()/p_0.sum()
            
            p0b = numpy.maximum(p_0 - 0.5, 0)
            z0 = (p0b*z_).sum()/p0b.sum()
            
            p_1 = 1.0*self.image.data[x,y,:,1].squeeze()
            p_1 -= p_1.min()
            p_1 /= p_1.max()
            
            p1b = numpy.maximum(p_1 - 0.5, 0)
            z1 = (p1b*z_).sum()/p1b.sum()
            
            dz = z1 - z0

            print(('z0: %f, z1: %f, dz: %f' % (z0,z1,dz)))
            
            pylab.plot(z_, p_0)
            pylab.plot(z_, p_1)
            pylab.vlines(z0, 0, 1)
            pylab.vlines(z1, 0, 1)
            pylab.figtext(.7,.7, 'dz = %3.2f' % dz)

def myplot_setlabel(xlabel=None,ylabel=None,title=None, label=None, xy=(0,0), ax=None, labsize=15,rightticks=False):
    import matplotlib as mpl
    mpl.rcParams['font.size'] = labsize+0.
#    mpl.rcParams['font.family'] = 'serif'# New Roman'
    mpl.rcParams['font.serif'] = 'Bitstream Vera Serif'

    mpl.rcParams['axes.labelsize'] = labsize+1.
    print labsize+1
    mpl.rcParams['xtick.labelsize'] = labsize+0.
    mpl.rcParams['ytick.labelsize'] = labsize+0.

    if label:
        print "######################## LABELS HERE##########################"
        
        if xy==(0,0):
            xy=(0.3,0.90)
        if not ax:
            print "WARNING: no axix, cannot place label"
            pl.figtext(xy[0],xy[1],label, fontsize=labsize)
        else:
            pl.text(xy[0],xy[1],label, transform=ax.transAxes, fontsize=labsize)
            if rightticks:
                ax.yaxis.tick_right()
    pl.xlabel(xlabel, fontsize=labsize+1)
    pl.ylabel(ylabel, fontsize=labsize+1)
    if title:
        pl.title(title)

    def dovis(self):
        """
        Do runtime visualization. 
        """

        pylab.clf()

        phi = self.cc_data.get_var("phi")

        myg = self.cc_data.grid

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

        pylab.xlabel("x")
        pylab.ylabel("y")
        pylab.title("phi")

        pylab.colorbar()
        
        pylab.figtext(0.05,0.0125, "t = %10.5f" % self.cc_data.t)

        pylab.draw()

    def buildAtmos(self, secz, xlim=[300, 1100], doPlot=False):
        """Generate the total atmospheric transmission profile at this airmass, using the coefficients C."""
        # Burke paper says atmosphere put together as 
        # Trans_total (alt/az/time) = Tgray * (e^-Z*tau_aerosol(alt/az/t)) * 
        #         * (1 - C_mol * BP(t)/BPo * A_mol(Z))  -- line 2
        #         * (1 - sqrt(C_mol * BP(t)/BPo) * A_mol(Z))  -- 3
        #         * (1 - C_O3 * A_O3(A) )
        #         * (1 - C_H2O(alt/az/time) * A_H2O(Z))
        # Tau_aerosol = trans['aerosol'] ... but replace with power law (because here all 1's)
        #  typical power law index is about tau ~ lambda^-1
        # A_mol = trans['O2']
                
        # secz = secz of this observation
        # wavelen / atmo_templates == building blocks of atmosphere, with seczlist / atmo_ind keys
        # C = coeffsdictionary = to, t1, t2, alpha0 (for aerosol), C_mol, BP, C_O3, C_H2O  values    
        if (abs(secz - self.secz) > interplimit):
            print "Generating interpolated atmospheric absorption profiles for this airmass %f" %(secz)
            self.interpolateSecz(secz)

        BP0 = 782 # mb
        # set aerosol appropriately with these coefficients
        self.atmo_abs['aerosol'] = 1.0 - numpy.exp(-secz * (self.C['t0'] + self.C['t1']*0.0 + self.C['t2']*0.0) 
                                                   * (self.wavelen/675.0)**self.C['alpha'])
        # set total transmission, with appropriate coefficients
        self.trans_total = numpy.ones(len(self.wavelen), dtype='float')
        self.trans_total = self.trans_total * (1.0 - self.C['mol'] * self.C['BP']/BP0 * self.atmo_abs['rayleigh'])  \
                      * ( 1 - numpy.sqrt(self.C['mol'] * self.C['BP']/BP0) * self.atmo_abs['O2']) \
                      * ( 1 - self.C['O3'] * self.atmo_abs['O3']) \
                      * ( 1 - self.C['H2O'] * self.atmo_abs['H2O']) \
                      * ( 1 - self.atmo_abs['aerosol'])
        # now we can plot the atmosphere
        if doPlot:
            pylab.figure()
            pylab.subplot(212)
            colorindex = 0
            for comp in self.atmo_ind:
                pylab.plot(self.wavelen, self.atmo_abs[comp], colors[colorindex], label='%s' %(comp))
                colorindex = self._next_color(colorindex)        
            leg =pylab.legend(loc=(0.88, 0.3), fancybox=True, numpoints=1, shadow=True)
            ltext = leg.get_texts()
            pylab.setp(ltext, fontsize='small')
            coefflabel = ""
            for comp in ('mol', 't0', 'alpha', 'O3', 'H2O'):
                coefflabel = coefflabel + "C[%s]:%.2f  " %(comp, self.C[comp])
                if (comp=='alpha') | (comp=='mol'):
                    coefflabel = coefflabel + "\n"
            pylab.figtext(0.2, 0.35, coefflabel, fontsize='small')
            pylab.xlim(xlim[0], xlim[1])
            pylab.ylim(0, 1.0)
            pylab.xlabel("Wavelength (nm)")
            pylab.subplot(211)
            pylab.plot(self.wavelen, self.atmo_trans[self.seczToString(1.2)]['comb'], 'r-', label='Standard X=1.2 (no aerosols)')
            pylab.plot(self.wavelen, self.trans_total, 'k-', label='Observed')
            leg = pylab.legend(loc=(0.12, 0.05), fancybox=True, numpoints=1, shadow=True)
            ltext = leg.get_texts()
            pylab.setp(ltext, fontsize='small')
            pylab.xlim(xlim[0], xlim[1])
            pylab.ylim(0, 1.0)
            pylab.title("Example Atmosphere at X=%.2f" %(secz))
        return

def plotData():
    marr = fetchData()
    textsize = 18

    yticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], xt, size=textsize)
    ylabel("Periode [s]", size=textsize)
    xticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], yt, size=textsize)
    xlabel("Ausstattung [%]", size=textsize)
    title(
        "Abweichung der Geschwindigkeit zwischen FCD und des simulierten Verkehrs", size=textsize)
    # title("Relative Anzahl erfasster Kanten", size=textsize)
    figtext(0.7865, 0.92, '[%]', size=textsize)

    # levels=arange(mmin-mmin*.1, mmax+mmax*.1, (mmax-mmin)/10.))
    contourf(marr, 50)
    # set fontsize and ticks for the colorbar:
    if showVal == EDGENO:
        cb = colorbar(ticks=[0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100])
    else:
        # cb = colorbar(ticks=range(17))
        cb = colorbar()
    for t in cb.ax.get_yticklabels():  # set colorbar fontsize of each tick
        t.set_fontsize(textsize)

    show()

	def feature_label(self, data):
		f = data[0]
		s = data[1]
		base_x = 0.88
		base_y = 0.74
		for i in range(0, len(data[0])):
			 pl.figtext(base_y, base_x-0.05*i, '{} -> {}'.format(f[i], s[i]))

 def gen_hrf(self):
     m_A, q_Z, mu_k, m_H, sigma_k, width, height, hrf0 = self.analy.Vbjde() 
     fgs = self.ConditionalNRLHist(m_A, q_Z)
     MMin = -1.0  # Y.min()
     MMax = 1.0  # Y.max()
     pas = (MMax - MMin) / 100
     xx = arange(MMin, MMax, pas)
     nf = 1
     g0 = self.gaussian(xx, mu_k[0][0], sigma_k[0][0])
     g1 = self.gaussian(xx, mu_k[0][1], sigma_k[0][1])
     print (g0, g1)
     fgs.insert(0, figure((self.nf + 1) * 123))
     title("Fonction de reponse", fontsize='xx-large')
     figtext(0.2, 0.04,
             'bande = ' + str(self.bande)+
             ' beta =' + str(self.beta) +
             ' sigma = ' + str(self.sigmaH) +
             ' pl = ' + str(self.pl) +
             ' dt = ' + str(self.dt) +
             ' thrf = ' + str(self.Thrf),
             #'mu_k = '+ str(self.mu_k) +
             #'sigma_k = '+ str(self.sigma_k),
            fontsize='x-large')
     plot(m_H)
     if self.shower == 1:
         show()
     return fgs

def main(args):
    if len(args) != 3:
        print "Usage: python plots.py <results file> <start time> <end time>"
        sys.exit()
    fh = open(args[0], "r")
    l = fh.readline().split()
    N = int(l[0])
    T = int(l[1])
    step = int(l[2])
    mean_c = float(l[3])
    ts = int(args[1]) / step
    td = int(args[2]) / step + step
    x = range(0, T + step, step)
    cl = []
    for line in fh.readlines():
        c = float(line)
        cl.append(c)
    fh.close()

    # Time evolution of the fraction of cooperators.
    pylab.figure(1, figsize = (7, 4.5), dpi = 500)
    pylab.xlabel(r"$t$")
    pylab.ylabel(r"$x$")
    pylab.plot(x[ts:td], cl[ts:td], "#000000", alpha = 0.6, linewidth = 2.0)
    pylab.figtext(0.82, 0.85, r"$x_\infty = %4.3f$" %(mean_c), 
                  ha = 'center', va = 'center', 
                  bbox = dict(facecolor = 'white', edgecolor = 'black'))
    pylab.xlim(int(args[1]), int(args[2]))
    pylab.ylim(0, 1)
    ax = pylab.gca()
    ax.xaxis.major.formatter.set_powerlimits((0,0))
    pylab.savefig("plot.pdf", format = "pdf")
    pylab.close(1)

	def add_label(self, feature_result, min_feature_label):
		base_x = 0.88
		base_y = 0.74
		for i in range(0, len(feature_result)):
			x1 = feature_result[i]
			x2 = min_feature_label[i]
			pl.figtext(base_y, base_x-0.05*i, '{} -> [{}]'.format(x1, x2))

 def showImage(self, xlim=None, ylim=None, clims=None, cmap=None, copy=False, stats=True):
     if copy:
         # useful if you're going to apply a hanning filter or something later, but
         # want to keep an imge of the original
         image = numpy.copy(self.image)
     else:
         image = self.image
     pylab.figure()
     pylab.title('Image')
     if xlim == None:
         x0 = 0
         x1 = self.nx
     else:
         x0 = xlim[0]
         x1 = xlim[1]
     if ylim == None:
         y0 = 0
         y1 = self.ny
     else:
         y0 = ylim[0]
         y1 = ylim[1]
     if clims == None:
         pylab.imshow(image, origin='lower', cmap=cmap)
     else:
         pylab.imshow(image, origin='lower', vmin=clims[0], vmax=clims[1], cmap=cmap)
     pylab.xlabel('X')
     pylab.ylabel('Y')
     cb = pylab.colorbar()
     clims = cb.get_clim()
     pylab.xlim(x0, x1)
     pylab.ylim(y0, y1)
     if stats:
         statstxt = 'Mean/Stdev/Min/Max:\n %.2f/%.2f/%.2f/%.2f' %(numpy.mean(image), numpy.std(image), image.min(), image.max())
         pylab.figtext(0.75, 0.03, statstxt)
     return clims

    def display_settings(self):
        from pylab import figtext
        assert self.search

        search = self.search
        figtext(.01,.99,search.str_graph_fixed(), va='top',fontsize='medium')
        figtext(.35, .99, search.str_graph_used(), va='top',fontsize='medium')

    def test_varying_inclination(self):
        #""" Test that the waveform is consistent for changes in inclination
        #"""
        sigmas = []
        incs = numpy.arange(0, 21, 1.0) * lal.PI / 10.0

        for inc in incs:
            # WARNING: This does not properly handle the case of SpinTaylor*
            # where the spin orientation is not relative to the inclination
            hp, hc = get_waveform(self.p, inclination=inc)
            s = sigma(hp, low_frequency_cutoff=self.p.f_lower)        
            sigmas.append(s)
         
        f = pylab.figure()
        pylab.axes([.1, .2, 0.8, 0.70])   
        pylab.plot(incs, sigmas)
        pylab.title("Vary %s inclination, $\\tilde{h}$+" % self.p.approximant)
        pylab.xlabel("Inclination (radians)")
        pylab.ylabel("sigma (flat PSD)")
        
        info = self.version_txt
        pylab.figtext(0.05, 0.05, info)
        
        if self.save_plots:
            pname = self.plot_dir + "/%s-vary-inclination.png" % self.p.approximant
            pylab.savefig(pname)

        if self.show_plots:
            pylab.show()
        else:
            pylab.close(f)

        self.assertAlmostEqual(sigmas[-1], sigmas[0], places=7)
        self.assertAlmostEqual(max(sigmas), sigmas[0], places=7)
        self.assertTrue(sigmas[0] > sigmas[5])

def do_panel(condition, refcondition, figtitle, cursor, band='r'):
    cmd = "select a.BT, f.flag, a.r_bulge,  a.n_bulge, a.ba_bulge, a.ba_disk, c.z from CAST as c, Flags_catalog as f, {band}_band_serexp as a, gz2_flags as z where a.galcount = c.galcount and a.galcount = f.galcount and a.galcount = z.galcount and f.band='{band}' and f.model='serexp' and f.ftype ='u' and {condition};"
#a.r_bulge,
    BT, flags, r_bulge, n_bulge,ba_bulge,ba_disk, zspec =cursor.get_data(cmd.format(condition=condition, band=band))
    ref_BT,ref_flag, ref_rb, ref_nb, ref_ba,ref_badisk, ref_zspec = cursor.get_data(cmd.format(condition=refcondition, band=band))
    BT = np.array(BT)
    flags = np.array(flags, dtype=int)
    r_bulge = np.array(r_bulge)
    n_bulge= np.array(n_bulge)
    ba_bulge= np.array(ba_bulge)
    ba_disk= np.array(ba_disk)
    zspec = np.array(zspec)

    ref_BT = np.array(ref_BT)
    ref_flag = np.array(ref_flag, dtype=int)
    ref_rb = np.array(ref_rb)
    ref_nb= np.array(ref_nb)
    ref_ba= np.array(ref_ba)
    ref_badisk= np.array(ref_badisk)
    ref_zspec = np.array(ref_zspec)


    panel_plot(BT, flags, r_bulge, n_bulge,ba_bulge,ba_disk, 
               ref_BT,ref_flag, ref_rb, ref_nb, ref_ba,ref_badisk, 
               zspec, ref_zspec)
    pl.figtext(0.5, 0.95, figtitle)
    print figtitle+' ',BT.size, ' objects' 
    pl.savefig('bar_params_serexp_%s.eps' %band)
    #pl.savefig('bar_z_serexp.eps')
    return

    def plot_Rnl(self,filename=None):
        """ Plot radial wave functions with matplotlib.
        
        filename:  output file name + extension (extension used in matplotlib)
        """
        if pl==None:
            raise AssertionError('pylab could not be imported')
        rmax = data[self.symbol]['R_cov']/0.529177*3
        ri = np.where( self.rgrid<rmax )[0][-1]
        states=len(self.list_states())
        p = np.ceil(np.sqrt(states)) #p**2>=states subplots
        
        fig=pl.figure()        
        i=1
        # as a function of grid points
        for n,l,nl in self.list_states():
            ax=pl.subplot(2*p,p,i)
            pl.plot(self.Rnlg[nl])
            pl.yticks([],[])
            pl.xticks(size=5)
            
            # annotate
            c = 'k'
            if nl in self.valence: 
                c='r'
            pl.text(0.5,0.4,r'$R_{%s}(r)$' %nl,transform=ax.transAxes,size=15,color=c)
            if ax.is_first_col():
                pl.ylabel(r'$R_{nl}(r)$',size=8)
            i+=1
            
        # as a function of radius
        i = p**2+1
        for n,l,nl in self.list_states():
            ax=pl.subplot(2*p,p,i)
            pl.plot(self.rgrid[:ri],self.Rnlg[nl][:ri])
            pl.yticks([],[])
            pl.xticks(size=5)
            if ax.is_last_row():
                pl.xlabel('r (Bohr)',size=8)

            c = 'k'
            if nl in self.valence: 
                c='r'
            pl.text(0.5,0.4,r'$R_{%s}(r)$' %nl,transform=ax.transAxes,size=15,color=c)
            if ax.is_first_col():
                pl.ylabel(r'$R_{nl}(r)$',size=8)
            i+=1
        

        file = '%s_KSAllElectron.pdf' %self.symbol
        #pl.rc('figure.subplot',wspace=0.0,hspace=0.0)
        fig.subplots_adjust(hspace=0.2,wspace=0.1)
        s=''
        if self.confinement!=None:
            s='(confined)'
        pl.figtext(0.4,0.95,r'$R_{nl}(r)$ for %s-%s %s' %(self.symbol,self.symbol,s))
        if filename is not None:
            file = filename
        pl.savefig(file)

 def putText( self, text ):
   """
   Puts some text into an otherwise empty plot.
   @param text: text to put in the empty plot
   """
   newText = ''
   for i in range( int(len(text)/60.0)+1):
     newText+=text[60*i:60*i+60]+'\n'
   figtext(0.15,0.15, newText)

def get_GB_cls_metrics(data_fh,info):
    """
    Get the metrics of Gradient Boost classification models

    :param data_classi_fh: path to file containing Classification training data
    """
    from pylab import figtext
    try:
        dpkl=read_pkl(data_fh)
    except:
        return False
    if not 'gs_cv' in dpkl.keys():
        return False
    dXy=dpkl['dXy_final']
    ycol=dpkl['ycol']
    gs_cv=dpkl['gs_cv']
    feat_imp = dpkl['feat_imp']

    Xcols=[c for c in dXy.columns.tolist() if c!=ycol]
    est=gs_cv.best_estimator_
    X=dXy.loc[:,Xcols].as_matrix()
    y=dXy.loc[:,ycol].as_matrix()        

    #partial dep 
    plot_type='partial_dep'
    plot_fh='%s/data_ml/%s.%s.pdf' % (info.prj_dh,plot_type,basename(data_fh))
    if not exists(plot_fh):
        feats_indi=[s for s in dpkl['feat_imp'].head(6).index.tolist() if not ((') ' in s) and (' (' in s))]
        features=[Xcols.index(f) for f in feats_indi]
        feature_names=linebreaker(Xcols)
        from sklearn.ensemble.partial_dependence import plot_partial_dependence
        fig, axs = plot_partial_dependence(est, X, features,#[[features[1],features[2]]],
                                           feature_names=feature_names,
                                           n_jobs=int(info.cores), grid_resolution=50,
                                           n_cols=2,
                                           line_kw={'color':'r'},
                                          figsize=[7,9])
        figtext(0.9,-0.2,'AUC = %.2f' % gs_cv.best_score_,ha='right',color='b')
        saveplot(plot_fh,form='pdf',tight_layout=False)
    
    #relimp
    plot_type='featimps'
    plot_fh='%s/data_ml/%s.%s.pdf' % (info.prj_dh,plot_type,basename(data_fh))
    if not exists(plot_fh):
        featst=10
        fig=plt.figure(figsize=(3,featst*0.75))
        # fig = plt.figure(figsize=(8,featst*0.25))#figsize=(11,5))
        ax=plt.subplot(111)
        feat_imp=feat_imp.sort_values(by='Feature importance',ascending=True)
        feat_imp.index=linebreaker(feat_imp.index, break_pt=30)
        feat_imp.tail(featst).plot(kind='barh',ax=ax, color='red')
        ax.set_xlabel('Feature Importance')
        ax.legend([])    
        figtext(0.9,-0.2,'AUC = %.2f' % gs_cv.best_score_,ha='right',color='b')
        saveplot(plot_fh,form='pdf',tight_layout=False)

    def test_varying_orbital_phase(self):
        #"""Check that the waveform is consistent under phase changes
        #"""
        
        if self.p.approximant in td_approximants():
            sample_attr = 'sample_times'
        else:
            sample_attr = 'sample_frequencies'   
            
        f = pylab.figure()
        pylab.axes([.1, .2, 0.8, 0.70])
        hp_ref, hc_ref = get_waveform(self.p, coa_phase=0)
        pylab.plot(getattr(hp_ref, sample_attr), hp_ref.real(), label="phiref")
       
        hp, hc = get_waveform(self.p, coa_phase=lal.PI/4)
        m, i = match(hp_ref, hp)
        self.assertAlmostEqual(1, m, places=2)
        o = overlap(hp_ref, hp)
        pylab.plot(getattr(hp, sample_attr), hp.real(), label="$phiref \pi/4$")
        
        hp, hc = get_waveform(self.p, coa_phase=lal.PI/2)
        m, i = match(hp_ref, hp)
        o = overlap(hp_ref, hp)
        self.assertAlmostEqual(1, m, places=7)
        self.assertAlmostEqual(-1, o, places=7)
        pylab.plot(getattr(hp, sample_attr), hp.real(), label="$phiref \pi/2$")
        
        hp, hc = get_waveform(self.p, coa_phase=lal.PI)
        m, i = match(hp_ref, hp)
        o = overlap(hp_ref, hp)
        self.assertAlmostEqual(1, m, places=7)
        self.assertAlmostEqual(1, o, places=7)
        pylab.plot(getattr(hp, sample_attr), hp.real(), label="$phiref \pi$")
        
        pylab.xlim(min(getattr(hp, sample_attr)), max(getattr(hp, sample_attr)))
        pylab.title("Vary %s oribital phiref, h+" % self.p.approximant)
        
        if self.p.approximant in td_approximants():
            pylab.xlabel("Time to coalescence (s)")
        else:
            pylab.xlabel("GW Frequency (Hz)") 

        pylab.ylabel("GW Strain (real part)")
        pylab.legend(loc="upper left")
        
        info = self.version_txt
        pylab.figtext(0.05, 0.05, info)
        
        if self.save_plots:
            pname = self.plot_dir + "/%s-vary-phase.png" % self.p.approximant
            pylab.savefig(pname)
        if self.show_plots:
            pylab.show()
        else:
            pylab.close(f)

def plot_peak_offsets(offsets, filebase, doshow=False):
    """Plots the results of the peak offsets calculated."""
    from pylab import figure, clf, xlim, ylim, savefig, plot, text, show, figtext, subplot, title
    #@type inst Instrument
    inst = instrument.inst
    numperpage = 6

    # Initialize some stats
    rms = 0
    rms_wl = 0
    #@type po PeakOffset
    for po in offsets:
        #Square of the error distance
        error = (po.measured[0]- po.predicted[0])**2 + (po.measured[1]- po.predicted[1])**2
        rms += error
        error = (po.wavelength_measured - po.wavelength_predicted)**2
        rms_wl += error
    # Now do the root-mean
    rms = (rms/ len(offsets))**0.5
    rms_wl = (rms_wl/ len(offsets))**0.5
    print "Peak offsets RMS error is ", rms
    print "Peak offsets RMS wavelength error is ", rms_wl

    #@type det FlatDetector
    for (det_num, det) in enumerate(inst.detectors):
        if det_num % numperpage == 0:
            figure(det_num/numperpage, figsize=[8, 10])
            clf()
            figtext(0.5, 0.95, "Offset (black line) from predicted peak positions (red dot); wavelength in angstroms", horizontalalignment='center', verticalalignment='top')
        ax = subplot(3, 2, det_num % numperpage+1)
        #Set the axes font sizes
        for xlabel_i in ax.get_xticklabels() + ax.get_yticklabels() :
            xlabel_i.set_fontsize(10)
        #@type po PeakOffset
        for po in offsets:
            if po.det_num == det_num:
                x = [po.measured[0], po.predicted[0]]
                y = [po.measured[1], po.predicted[1]]
                plot(po.predicted[0], po.predicted[1], 'r.')
                plot(x, y, '-k')
                text(po.predicted[0], po.predicted[1], ' %.1f' % po.wavelength_measured, size=5, verticalalignment='center')
        xlim( -det.width/2, det.width/2)
        ylim( -det.height/2, det.height/2)
        #axis('equal')
        title('Detector %s' % det.name)
    #-- Save to files --
    for i in xrange((len(inst.detectors) + numperpage-1) / numperpage):
        figure(i)
        savefig( filebase + "_%d.pdf" % i, papertype="letter")
    #-- combine --
    os.system("pdftk %s_*.pdf cat output %s.pdf" % (filebase, filebase))

    if doshow:
        show()

def showdp(fsz=16):
    """Utility function to set default parameters for DOS plots."""
    pl.xlabel("E-E$_\mathrm{f}$ (eV)", size=fsz)
    pl.figtext(0.03, 0.45, "LDOS", rotation="vertical", size=fsz)
    loc, lab = pl.xticks()
    lab.set_size = fsz
    loc, lab = pl.yticks()
    lab.set_size = fsz
    pl.legend()
    pl.subplots_adjust(hspace=0.0)
    pl.show()