def plot_data_objs(data_to_plot_list, figure_name, same=False):
    rows = 2
    cols = 2
    topo = figure_name.split(':')[0]
    # print "topo:",topo
    # print "figname:",figure_name
    if same:
        fig_id = figure_name.split(':')[-1]
        fig, ax = plt.subplots(rows, cols, num=fig_id)
    else:
        fig, ax = plt.subplots(rows, cols)
    for i in range(rows):
        for j in range(cols):
            if len(data_to_plot_list) > i*rows+j:
                obj = data_to_plot_list[i*rows+j]
                if obj.scale == 'log' and max(obj.data) <= 0:
                    ydata = map(lambda x: -x, obj.data)
                else:
                    ydata = obj.data
                ax[i][j].plot(range(len(obj.data)), ydata, label=topo, linewidth=2.0)
                ax[i][j].set_yscale(obj.scale)
                ax[i][j].set_title(obj.label)
                if same:
                    ax[i][j].legend(loc='lower left', shadow=True)
                    ax[i][j].legend().set_visible(False)
                    # fig.legend(loc='lower left', shadow=True)
    if fig._suptitle is None:
        suptitle = ':'.join(figure_name.split(':')[1:])
        fig.suptitle(suptitle, fontsize=18)

	def plotStrainHeatmap(self, alldata):
		img = np.zeros((len(alldata[0, :, 0]), len(alldata[0, 0, :]), 4),  dtype=np.uint8)

		# for i in range(len(alldata[:, 0, 0])/2+1):
		# 	rawdata = alldata[i, :, :]
		# 	print i, np.max(rawdata)
		# 	img[:, :, 3] = 255  # -20*i
		# 	img[:, :, 1] = (255-i*10)*rawdata/np.max(rawdata)
		# 	img[:, :, 2] = 0  # 255*rawdata/np.max(rawdata)
		# 	img[:, :, 0] = 0  # 255*rawdata/np.max(rawdata)
		#
		# for i in range(len(alldata[:, 0, 0])/2):
		# 	print i+len(alldata[:, 0, 0])/2+1, np.max(rawdata)
		# 	rawdata = alldata[i+len(alldata[:, 0, 0])/2+1, :, :]
		# 	img[:, :, 3] = 255  # -20*i
		# 	img[:, :, 1] = 0
		# 	img[:, :, 2] = (180+i*10)*rawdata/np.max(rawdata)
		# 	img[:, :, 0] = (180+i*10)*rawdata/np.max(rawdata)

		img = alldata[10, :, :]

		print np.mean(img), np.max(img)
		print np.mean(alldata[2, :, :]), np.max(alldata[2, :, :])

		linestart = [999, 0]
		linestop = [0, 999]

		z, length = self.getProjection(img, linestart[0], linestart[1], linestop[0], linestop[1])

		fig, ax = plt.subplots()
		fig2, ax2 = plt.subplots()

		ax2.plot(z)

		ax.imshow(img)

 def plotVolume(self):
     full = self.revolution(1000)
     volu = self.volume(full)
     
     import matplotlib as mpl
     #mpl.use('Qt4Agg')
     from matplotlib.pyplot import plot, show
     import matplotlib.pylab as plt
     
     if self.DEBUG:
         fig, axs = plt.subplots(2, 1, sharex=True)
         rev = self.revolution(2000)
         pos  = self._position(rev)
         axs[0].plot(rev*180/pi,pos*100)
         axs[0].plot(self.TDC()*180/pi,self._position(self.TDC())*100,'o')
         iMin  = np.where(pos==pos.min())
         axs[0].plot(rev[iMin]*180/pi,self._position(rev[iMin])*100,'o')
         iMax  = np.where(pos==pos.max())
         axs[0].plot(rev[iMax]*180/pi,self._position(rev[iMax])*100,'o')
         axs[0].set_ylabel(r'Piston position (cm)')
         ax = axs[1]
         self.autolog("Position: ", str(pos.min()), str(pos.max()), str(self.stroke()))
         self.autolog("Volume:    ",str(volu.min()), str(volu.max()))
     else:
         fig, ax = plt.subplots(1, 1)
         
     ax.plot(full*180/pi,volu*1e6)
     ax.plot(self.TDC()*180/pi,self.volume(self.TDC())*1e6,'o')
     iMin  = np.where(volu==volu.min())
     ax.plot(full[iMin]*180/pi,self.volume(full[iMin])*1e6,'o')
     iMax  = np.where(volu==volu.max())
     ax.plot(full[iMax]*180/pi,self.volume(full[iMax])*1e6,'o')
     ax.set_xlabel(r'Crankshaft angle $\theta$ (deg)')
     ax.set_ylabel(r'Cylinder volume $V$ (cm$^3$)')
     show()

	def plotStrain(self, strainpic, imgarray, label):
		# import matplotlib.ticker as ticker
		# sns.set_context("talk")

		print "strainpic dimensions: " + str(np.shape(strainpic))

		# gradient = self.adjustGradient()

		figstrain, axstrain = plt.subplots(2, 1)

		strainpic_adjusted = strainpic - np.mean(strainpic)
		strainpic_adjusted[strainpic_adjusted > 0.00004] = 0.00004
		strainpic_adjusted[strainpic_adjusted < -0.00004] = -0.00004

		im = axstrain[0].imshow(strainpic_adjusted, cmap="BrBG")
		# im = axstrain[0].imshow(strainpic+gradient, cmap="BrBG")
		im2 = axstrain[1].imshow(imgarray[len(imgarray[:, 0, 0])/2-4, :, :], cmap="Greens")

		axstrain[1].set_title("%g %g %g %g" % (self.roi[0], self.roi[1], self.roi[2], self.roi[3]))
		axstrain[0].set_title(r'$\epsilon_{220}$')

		def fmt(x,  pos):
			a,  b = '{:.2e}'.format(x).split('e')
			b = int(b)
			return r'${} \times 10^{{{}}}$'.format(a,  b)

		figstrain.subplots_adjust(right=0.8)
		cbar_ax1 = figstrain.add_axes([0.85,  0.55,  0.02,  0.35])
		cbar_ax2 = figstrain.add_axes([0.85,  0.1,  0.02,  0.35])
		clb = figstrain.colorbar(im, cax=cbar_ax1)  # , format=ticker.FuncFormatter(fmt))
		figstrain.colorbar(im2, cax=cbar_ax2)

		linestart = [100, 50]
		linestop = [100, 350]
		clb.set_clim(-0.00004, 0.00004)
		axstrain[0].autoscale(False)
		axstrain[0].plot([linestart[0], linestop[0]], [linestart[1], linestop[1]])

		z, length = self.getProjection(strainpic, linestart[0], linestart[1], linestop[0], linestop[1], 500)
		f3, ax3 = plt.subplots()
		linerange = np.linspace(0, 90*length/1000, len(z))
		ax3.plot(linerange, z)
		ax3.set_ylabel(r'Strain [$\Delta\theta/\theta$]')
		ax3.set_xlabel(r'[$\mu m$]')

		# np.save(self.directory + '/strainmap_array.txt', strainpic)
		f3.savefig(self.directory + '/strainmap_line.pdf')
		figstrain.savefig(self.directory + '/strainmap_%s.pdf' % str(label))
		# f4, ax4 = plt.subplots()
		# strain = np.reshape(strainpic, len(strainpic[:, 0])*len(strainpic[0, :]))
		# # strain[strain<-0.0005] = 0
		# # strain[strain>0.0005] = 0
		# # sns.distplot(strain, kde=False,  rug=False)
		# ax4.set_xlim(np.min(strain)-abs(0.1*np.min(strain)), np.max(strain)+0.1*np.max(strain))
		# # ax4.set_xlim(-0.0004,0.0004)
		# ax4.set_xlabel(r'$\theta$ offset [$^o$]')
		# ax4.set_title('Strain distribution')
		# f4.savefig(self.directory + '/straindistribution.pdf')
		return figstrain, axstrain

def plot_norm(samples, fields=None, filename=None):
    # get time from timestamp and sample time
    t = samples.bicycle.dt.mean() * samples.ts
    n = t.shape[0]

    if fields is None:
        # Set default to be fields that are not scalar with data available for
        # at least 10% of the timerange
        fields = []
        for name in samples.dtype.names:
            mask = samples.mask[name]
            if len(mask.shape) < 2:
                continue
            if reduce(mul, mask.shape[1:], 1) == 1:
                continue
            if np.count_nonzero(mask) < n/10:
                fields.append(name)
    if isinstance(fields, str):
        fields = (fields,)
    n = len(fields)
    if n > 6:
        color = sns.husl_palette(n)
    else:
        color = sns.color_palette('muted', n)

    if n > 1:
        fig, axes = plt.subplots(math.ceil(n/2), 2)
        axes = np.ravel(axes)
        if len(axes) > n:
            axes[-1].axis('off')
    else:
        fig, ax = plt.subplots()
        axes = [ax]

    for n, f in enumerate(fields):
        ax = axes[n]
        X = samples.__getattribute__(f)
        x = np.linalg.norm(X, axis=(1, 2))

        ax.set_xlabel('{} [{}]'.format('time', unit('time')))
        ax.plot(t, x, color=color[n], label=f)
        ax.legend()

    if n > 0:
        title = 'Norms'
    else:
        ax.legend().remove()
        field_parts = fields[0].split('.')
        title = ' '.join([f.title() for f in field_parts[:-1]] +
                         field_parts[-1:])
        title = 'Norm of ' + title
        axes = ax
    _set_suptitle(fig, title, filename)
    return fig, axes

def apply_unrot(filename):

    import read_idb as ri
    import dbutil as db
    import copy
    from util import lobe, Time
    import matplotlib.pylab as plt
    import numpy as np
    blah = np.load('/common/tmp/Feed_rotation/20170702121949_delay_phase.npz')
    dph = blah['dph']
    fghz = blah['fghz']
    out = ri.read_npz([filename])
    nbl, npol, nfrq, nt = out['x'].shape
    # Correct data for phase
    #n = [0,0,0,1,1,0,1,0,1,1,0,0,0]
    for i in range(13):
        a1 = lobe(dph[i] - dph[13])
        a2 = -dph[13] + np.pi/2
        a3 = dph[i] - np.pi/2
        for j in range(nt):
            out['x'][ri.bl2ord[i,13],1,:,j] *= np.exp(1j*a1)
            out['x'][ri.bl2ord[i,13],2,:,j] *= np.exp(1j*a2) 
            out['x'][ri.bl2ord[i,13],3,:,j] *= np.exp(1j*a3)
        
    trange = Time(out['time'][[0,-1]],format='jd')
    times, chi = db.get_chi(trange)
    nskip = len(times)/nt
    chi = np.transpose(chi[::nskip+1])
    chi[[8,9,10,12]] = 0.0
    outp = copy.deepcopy(out)
    for i in range(nt):
        for k in range(13):
            outp['x'][ri.bl2ord[k,13],0] = out['x'][ri.bl2ord[k,13],0]*np.cos(chi[k,i]) + out['x'][ri.bl2ord[k,13],3]*np.sin(chi[k,i])
            outp['x'][ri.bl2ord[k,13],2] = out['x'][ri.bl2ord[k,13],2]*np.cos(chi[k,i]) + out['x'][ri.bl2ord[k,13],1]*np.sin(chi[k,i])
            outp['x'][ri.bl2ord[k,13],3] = out['x'][ri.bl2ord[k,13],3]*np.cos(chi[k,i]) - out['x'][ri.bl2ord[k,13],0]*np.sin(chi[k,i])
            outp['x'][ri.bl2ord[k,13],1] = out['x'][ri.bl2ord[k,13],1]*np.cos(chi[k,i]) - out['x'][ri.bl2ord[k,13],2]*np.sin(chi[k,i])
    amp0 = np.abs(np.sum(out['x'][ri.bl2ord[:,13]],3))
    amp2 = np.abs(np.sum(outp['x'][ri.bl2ord[:,13]],3))
    f, ax = plt.subplots(4,13)
    for i in range(13):
        for j in range(4):
            ax[j,i].cla()
            ax[j,i].plot(fghz, amp0[i,j],'.',color='lightgreen')
            ax[j,i].plot(fghz, amp2[i,j],'k.')
    ph0 = np.angle(np.sum(out['x'][ri.bl2ord[:,13]],3))
    ph2 = np.angle(np.sum(outp['x'][ri.bl2ord[:,13]],3))
    f, ax = plt.subplots(4,13)
    for i in range(13):
        for j in range(4):
            ax[j,i].cla()
            ax[j,i].plot(fghz, ph0[i,j],'.',color='lightgreen')
            ax[j,i].plot(fghz, ph2[i,j],'k.')
            

def makePlots(sampletitle, amplpic, midppic, fwhmpic):
    fig0, ax0 = plt.subplots(1, 1, dpi=150, figsize=[5, 5])
    # plt.tight_layout()
    fig1, ax1 = plt.subplots(1, 1, dpi=150, figsize=[5, 5])
    # plt.tight_layout()
    fig2, ax2 = plt.subplots(1, 1, dpi=150, figsize=[5, 5])
    # plt.tight_layout()

    # fig0.set_figsize_inches(7,7)
    # fig1.set_size_inches(7,7)
    # fig2.set_size_inches(7,7)

    im0 = ax0.imshow(amplpic[3:-3, 3:-3], cmap='jet', interpolation='None')
    im1 = ax1.imshow(fwhmpic[3:-3, 3:-3], cmap='BrBG', interpolation='None')
    im2 = ax2.imshow(midppic[3:-3, 3:-3] , cmap='BrBG', interpolation='None')

    ax0.set_title('AMPL')
    ax1.set_title('FWHM')
    ax2.set_title('MIDP')

    def fmt(x,  pos):
        a,  b = '{:.2e}'.format(x).split('e')
        b = int(b)
        return r'${} \times 10^{{{}}}$'.format(a,  b)

    fig0.subplots_adjust(right=0.8)
    fig1.subplots_adjust(right=0.8)
    fig2.subplots_adjust(right=0.8)

    cbar_ax0 = fig0.add_axes([0.85,  0.1,  0.05,  0.8])
    cbar_ax1 = fig1.add_axes([0.85,  0.1,  0.05,  0.8])
    cbar_ax2 = fig2.add_axes([0.85,  0.1,  0.05,  0.8])

    clb0 = fig0.colorbar(im0, cax=cbar_ax0)  # , format=ticker.FuncFormatter(fmt))
    clb1 = fig1.colorbar(im1, cax=cbar_ax1)  # , format=ticker.FuncFormatter(fmt))
    clb2 = fig2.colorbar(im2, cax=cbar_ax2)  # , format=ticker.FuncFormatter(fmt))

    # # clb0.set_clim(0., 200.)
    # clb2.set_clim(-0.1, 0.1)
    # clb1.set_clim(-0.03, 0.03)

    fig0.savefig('plots/%s-ampl.pdf' % (sampletitle))
    fig1.savefig('plots/%s-fwhm.pdf' % (sampletitle))
    fig2.savefig('plots/%s-midp.pdf' % (sampletitle))

    fig0.clf()
    fig1.clf()
    fig2.clf()

def exampleNetworks():
    names = [u'Yoachim, P', u'Bellm, E', u'Williams, B', u'Williams, B', u'Capelo, P']
    # add some caching so it only querries once.
    if not hasattr(exampleNetworks,'results'):
        exampleNetworks.results = [None for name in names]
        exampleNetworks.graphs = [None for name in names]

    years = [2007, 2011, 2002, 2010, 2012]
    texts = ['(a)', '(b)','(c)', '(d)','(e)']
    count = 1
    figs = []
    filenames = []
    for name,year,txt in zip(names,years,texts):
        fig,ax = plt.subplots()
        figDummy, axDummy = plt.subplots()
        phdA =  list(ads.SearchQuery(q=u'bibstem:*PhDT', author=name, year=year,
                                     database='astronomy'))[-1]
        if exampleNetworks.results[count-1] is None:
            result, graph = phdArticle2row(phdA, checkUSA=False, verbose=True, returnNetwork=True)
            exampleNetworks.results[count-1] = result
            exampleNetworks.graphs[count-1] = graph
        else:
            result = exampleNetworks.results[count-1]
            graph = exampleNetworks.graphs[count-1]
        years = []
        for node in graph.nodes():
            years.append(float(node[0:4]))

        years = np.array(years)
        # Make the graph repeatable
        pos = {}
        for i, node in enumerate(graph.nodes()):
            pos[node] = (years[i],i**2)
        layout = nx.spring_layout(graph, pos=pos)
        nx.draw_networkx(graph, pos=layout, ax=ax, node_size=100,
                         node_color=years, alpha=0.5, with_labels=False)
        #nx.draw_spring(graph, ax=ax, node_size=100,
        #               node_color=years, alpha=0.5, with_labels=False)
        mappableDummy = axDummy.scatter(years,years,c=years)
        cbar = plt.colorbar(mappableDummy, ax=ax, format='%i')
        cbar.set_label('Year')
        ax.text(.1,.8, txt, fontsize=24, transform=ax.transAxes)
        ax.set_axis_off()
        figs.append(fig)
        filenames.append('example_network_%i' %count)
        count += 1
        print result
    return figs, filenames

    def plot(self,frequency=True,phase=False,dB=True,cal=True,fig=[],ax=[],color='k'):
        """

        """
        if fig==[]:
            fig,ax=plt.subplots(8,self.Nt,sharex=True,sharey=True)
        if cal:
            H = self.Hcal
        else:
            H = self.H
        for iR in range(self.Nr):
            for iT in range(self.Nt):
                k = iR*4+iT
                if frequency:
                    if not phase:
                        if dB:
                            #ax[iR,iT].plot(H.x,20*np.log10(abs(H.y[k,:])),color=color)
                            ax[iR,iT].plot(H.x,20*np.log10(abs(H.y[iR,iT,:])),color=color)
                        else:
                            #ax[iR,iT].plot(H.x,abs(H.y[k,:]),color='k')
                            ax[iR,iT].plot(H.x,abs(H.y[iR,iT,:]),color='k')
                    else:
                        #ax[iR,iT].plot(H.x,np.unwrap(np.angle(H.y[k,:])),color=color)
                        ax[iR,iT].plot(H.x,np.unwrap(np.angle(H.y[iR,iT,:])),color=color)
                else:
                        ax[iR,iT].plot(self.h.x,abs(self.h.y[iR,iT,:]),color=color)
                if (iR==7):
                    ax[iR,iT].set_xlabel('f (GHz)')
                ax[iR,iT].set_title(str(iR+1)+'x'+str(iT+1))
        return(fig,ax)

def plotGetRetangle():
    """ Area selection from selected pen.
    """
    selRect = []
    if len(ds.EpmDatasetAnalysisPens.SelectedPens) != 1:
        sr.msgBox('EPM Python Plugin - Demo Tools', 'Please select a single pen before applying this function!', 'Warning')
        return 0
    epmData = ds.EpmDatasetAnalysisPens.SelectedPens[0].values
    y = epmData['Value'].copy()
    x = np.arange(len(y))
    fig, current_ax = pl.subplots()
    pl.plot(x, y, lw=2, c='g', alpha=.3)

    def line_select_callback(eclick, erelease):
        'eclick and erelease are the press and release events'
        x1, y1 = eclick.xdata, eclick.ydata
        x2, y2 = erelease.xdata, erelease.ydata
        print ("\n(%3.2f, %3.2f) --> (%3.2f, %3.2f)" % (x1, y1, x2, y2))
        selRect.append((int(x1), y1, int(x2), y2))

    def toggle_selector(event):
        if event.key in ['Q', 'q'] and toggle_selector.RS.active:
            toggle_selector.RS.set_active(False)
        if event.key in ['A', 'a'] and not toggle_selector.RS.active:
            toggle_selector.RS.set_active(True)
    toggle_selector.RS = RectangleSelector(current_ax, line_select_callback, drawtype='box', useblit=True, button=[1,3], minspanx=5, minspany=5, spancoords='pixels')
    pl.connect('key_press_event', toggle_selector)
    pl.show()
    return selRect

def plot_MheF(isotracks=None, labels=None, colors=None):
    """ plot the minimum initial mass for He Fusion """
    if isotracks is None:
        isotracks = ['isotrack/parsec/CAF09_MC_S13v3_OV0.3.dat',
                     'isotrack/parsec/CAF09_MC_S13v3_OV0.4.dat',
                     'isotrack/parsec/CAF09_MC_S13v3_OV0.5.dat',
                     'isotrack/parsec/CAF09_MC_S13v3_OV0.6.dat',
                     'isotrack/parsec/CAF09_S12D_NS_1TP.dat']
        isotracks = [os.path.join(os.environ['TRILEGAL_ROOT'], i)
                     for i in isotracks]

    if labels is None:
        labels = ['$\Lambda_c=0.3$',
                  '$\Lambda_c=0.4$',
                  '$\Lambda_c=0.5$',
                  '$\Lambda_c=0.6$',
                  '$S12D\_NS\_1TP$']
    if colors is None:
        colors = ['darkred', 'orange', 'navy', 'purple', 'k']

    fig, ax = plt.subplots()
    for i, isotrack in enumerate(isotracks):
        isot = trilegal.IsoTrack(isotrack)
        ax.plot(isot.Z, isot.mhefs, lw=2, label=labels[i], color=colors[i])
        ax.plot(isot.Z, isot.mhefs, 'o', color=colors[i])

    ax.grid()
    ax.set_xlim(0.001, 0.0085)
    ax.set_ylim(1.55, 2.05)
    return ax

def test_params():
    #x = np.linspace(.8, 1.2, 1e2)
    x = np.linspace(-.2, .2, 1e2)

    num = 5
    range_a = np.linspace(1, 2, num)
    range_b = np.linspace(1., 1.1, num)
    range_p = np.linspace(.1, .4, num)
    range_q = np.linspace(.1, .4, num)
    range_T = np.linspace(30, 365, num) / 365

    args_def = {'a' : range_a.mean(), 'b' : range_b.mean(),
                'p' : range_p.mean(), 'q' : range_q.mean(),
                'T' : range_T.mean()}

    ranges = {'a' : range_a, 'b' : range_b,
              'p' : range_p, 'q' : range_q, 'T' : range_T}

    fig, axes = plt.subplots(nrows = len(ranges), figsize = (6,12))
    for name, a in zip(sorted(ranges.keys()), axes):
        args = args_def.copy()
        for pi in ranges[name]:
            args[name] = pi
            f = GB2(**args).density(x)
            a.plot(x, f, label = pi)
        a.legend(title = name)
    plt.show()

def plot_integrated_colors(filenames, labels='Z'):
    if type(filenames) is str:
        filenames = [filenames]
        ax = None
        cols = ['k']
    else:
        fig, ax = plt.subplots()
        cols = brewer2mpl.get_map('Spectral', 'Diverging',
                                  len(filenames)).mpl_colors

    if labels == 'Z':
        fmt = '$Z=%.4f$'
        labels = [fmt % float(l.replace('.dat', '').split('Z')[1])
                  for l in filenames]
    else:
        print 'need to fix labels'
        labels = [''] * len(filenames)
    for i, filename in enumerate(filenames):
        data = rsp.fileIO.readfile(filename)
        ycol = 'V-K'
        xcol = 'Age'
        ax = rg.color_color(data, xcol, ycol, xscale='log', ax=ax,
                            plt_kw={'lw': 2, 'color': cols[i],
                                    'label': labels[i]})

    plot_cluster_data(ax)
    ax.legend(frameon=False, loc=0, numpoints=1)
    ax.set_xlabel(r'${\rm %s}$' % xcol, fontsize=20)
    ax.set_ylabel(r'${\rm %s}$' % ycol, fontsize=20)
    plt.tick_params(labelsize=16)
    return ax

def interev_mag(times, mags):
    r"""Function to plot interevent times against magnitude for given times
    and magnitudes.

    :type times: list of datetime
    :param times: list of the detection times, must be sorted the same as mags
    :type mags: list of float
    :param mags: list of magnitudes
    """
    l = [(times[i], mags[i]) for i in xrange(len(times))]
    l.sort(key=lambda tup: tup[0])
    times = [x[0] for x in l]
    mags = [x[1] for x in l]
    # Make two subplots next to each other of time before and time after
    fig, axes = plt.subplots(1, 2, sharey=True)
    axes = axes.ravel()
    pre_times = []
    post_times = []
    for i in range(len(times)):
        if i > 0:
            pre_times.append((times[i] - times[i - 1]) / 60)
        if i < len(times) - 1:
            post_times.append((times[i + 1] - times[i]) / 60)
    axes[0].scatter(pre_times, mags[1:])
    axes[0].set_title('Pre-event times')
    axes[0].set_ylabel('Magnitude')
    axes[0].set_xlabel('Time (Minutes)')
    plt.setp(axes[0].xaxis.get_majorticklabels(), rotation=30)
    axes[1].scatter(pre_times, mags[:-1])
    axes[1].set_title('Post-event times')
    axes[1].set_xlabel('Time (Minutes)')
    plt.setp(axes[1].xaxis.get_majorticklabels(), rotation=30)
    plt.show()

def plot_locality_regression(snps,cob,gene_limit=10):
    # Get degree and bootstrap degree
    log('Fetching Empirical Degree')
    degree = cob.locality(cob.refgen.candidate_genes(snps,gene_limit=gene_limit,chain=True)).sort('local')
    log('Fetching BS Degree')
    #bsdegree = pd.concat([cob.locality(cob.refgen.bootstrap_candidate_genes(snps,gene_limit=gene_limit,chain=True)) for x in range(50)]).sort('local')
    # get OLS for the bootstrapped degree 
    log('Fitting models')
    model = sm.OLS(degree['global'],degree.local)
    res = model.fit()
    std, iv_l, iv_u = wls_prediction_std(res)
    # plot the bootstrapped data
    fig,ax = pylab.subplots(figsize=(8,6)) 
    fig.hold(True)
    ax.set_xlim(0,max(degree.local))
    ax.set_ylim(0,max(degree['global']))
    # plot the bootstraps std
    # plot the true data
    log('Plotting Empirical')
    ax.plot(degree.local,degree['global'],'o',label='Empirical')
    log('Plotting Residuals')
    ax.plot(degree.local,res.fittedvalues,'--')
    ax.plot(degree.local,res.fittedvalues+2.5*std,'r--')
    ax.plot(degree.local,res.fittedvalues-2.5*std,'r--')
    ax.set_xlabel('Number Local Interactions')
    ax.set_ylabel('Number Global Interactions')
    log('Saving Figure')
    fig.savefig('{}_locality.png'.format(cob.name))

def avg_scores_plot(fit_results, predictor_variable):
    n_array = []
    r2_array = []
    rmse_array = []
    for key in fit_results:
        value = fit_results[key]
        i = 0
        sum_r2 = 0
        sum_rmse = 0
        while i < len(value):
            sum_r2 = float(sum_r2) + value[i][1]
            sum_rmse = float(sum_rmse) + value[i][2]
            i = i + 1
        avg_r2 = sum_r2/len(value)
        avg_rmse = sum_rmse/len(value)
        n_array.append(key)
        r2_array.append(avg_r2)
        rmse_array.append(avg_rmse)
        print 'For n = ' + str(key) + ': Average R^2 = ' + str(avg_r2) + ', Average RMSE = ' + str(avg_rmse)  
    print 'Minimum Average RMSE is '+str(min(rmse_array))
    #Plot Average Values to determine best fit
    f, ax = plt.subplots(2, sharex=True)
    ax[0].scatter(n_array, r2_array, color='r', label='Average R^2 Value')
    ax[0].set_ylabel('Average R^2 Value')
    ax[0].xaxis.grid()
    ax[0].yaxis.grid()
    ax[1].scatter(n_array, rmse_array, color='blue', label='Average RMSE Value')
    ax[1].set_ylabel('Average RMSE Value')
    ax[1].xaxis.grid()
    ax[1].yaxis.grid()
    plt.xlabel('Polynomial Fit Order')
    ax[0].set_title('Avg. R^2 & RMSE in ' + str(len(value))+ '-fold Cross Validation of ' + str(predictor_variable) + ' to No. 311 Incidents')
    plt.show()

def plot_dataframe_meshgrid(df, xaxis = 0, ax = None):
    axes_list = [df.index, df.columns]
    x_index = axes_list[xaxis]

    if xaxis !=0:
        df = df.swapaxes(0,1)
        y_index = axes_list[0]
    else:
        y_index = axes_list[1]

    z = df.values.transpose()
    x = np.repeat(np.array([x_index]), y_index.shape[0], axis = 0)
    y = np.repeat(np.array([y_index]), x_index.shape[0], axis = 0).transpose()

    if ax:
        a = ax
        f = a.get_figure()
    else:
        f,a = plt.subplots()

    pc = a.pcolormesh(x, y , z)


    if 'datetime' in df.index.dtype_str:
        f.autofmt_xdate()
    cb = f.colorbar(pc)
    a.set_xlabel(df.index.name)
    a.set_ylabel(df.columns.name)

    # nans, screw up the scaling, therefore ...
    if np.any(np.isnan(df.values)):
        values = df.values
        values = values[~ np.isnan(values)]
        pc.set_clim((values.min(),values.max()))
    return f,a,pc,cb

def make_lf(ghist, gbins, shist, sbins, gerrs=None, serrs=None,
            colors=['black', 'darkred'], label='', ax=None):

    if ax is None:
        fig, ax = plt.subplots(figsize=(8, 8))
    if not None in [gerrs, serrs]:
        plt_kw = {'drawstyle': 'steps-mid', 'color': colors[0],
                  'lw': 2}
        ax.errorbar(gbins[1:], ghist, yerr=gerrs, **plt_kw)
        plt_kw['color'] = colors[1]
        plt_kw['label'] = label
        plt_kw['alpha'] = 0.3
        ax.errorbar(sbins[1:], shist, yerr=serrs, **plt_kw)
    else:
        plt_kw = {'linestyle': 'steps-mid', 'color': colors[0],
                  'lw': 2}
        ax.plot(gbins[1:], ghist, **plt_kw)
        plt_kw['color'] = colors[1]
        plt_kw['label'] = label
        plt_kw['alpha'] = 0.3
        ax.plot(sbins[1:], shist, **plt_kw)

    ax.set_yscale('log')
    ax.tick_params(labelsize=16)

    ax.set_ylabel(r'${\rm Number\ of\ Stars}$', fontsize=20)
    return ax

def SVD_plot(SVStreams, SValues, stachans, title=False):
    r"""Function to plot the singular vectors from the clustering routines, one\
    plot for each stachan

    :type SVStreams: list of :class:Obspy.Stream
    :param SVStreams: See clustering.SVD_2_Stream - will assume these are\
            ordered by power, e.g. first singular vector in the first stream
    :type SValues: list of float
    :param SValues: List of the singular values corresponding to the SVStreams
    :type stachans: list
    :param stachans: List of station.channel
    """
    for stachan in stachans:
        print(stachan)
        plot_traces = [SVStream.select(station=stachan.split('.')[0],
                                       channel=stachan.split('.')[1])[0]
                       for SVStream in SVStreams]
        fig, axes = plt.subplots(len(plot_traces), 1, sharex=True)
        axes = axes.ravel()
        for i, tr in enumerate(plot_traces):
            y = tr.data
            x = np.linspace(0, len(y) * tr.stats.delta, len(y))
            axes[i].plot(x, y, 'k', linewidth=1.1)
            ylab = 'SV '+str(i+1)+'='+str(round(SValues[i] / len(SValues), 2))
            axes[i].set_ylabel(ylab, rotation=0)
            axes[i].yaxis.set_ticks([])
            print(i)
        axes[-1].set_xlabel('Time (s)')
        plt.subplots_adjust(hspace=0)
        if title:
            axes[0].set_title(title)
        else:
            axes[0].set_title(stachan)
        plt.show()
    return

    def check_models(self):
        temp = np.logspace(0, np.log10(600))
        num = len(self.available_models())

        fig, ax = plt.subplots(1)
        self.plotting_colours(num, fig, ax, repeats=2)

        for author in self.available_models():
            Nc, Nv = self.update(temp=temp, author=author)
            # print Nc.shape, Nv.shape, temp.shape
            ax.plot(temp, Nc, '--')
            ax.plot(temp, Nv, '.', label=author)

        ax.loglog()
        leg1 = ax.legend(loc=0, title='colour legend')

        Nc, = ax.plot(np.inf, np.inf, 'k--', label='Nc')
        Nv, = ax.plot(np.inf, np.inf, 'k.', label='Nv')

        plt.legend([Nc, Nv], ['Nc', 'Nv'], loc=4, title='Line legend')
        plt.gca().add_artist(leg1)

        ax.set_xlabel('Temperature (K)')
        ax.set_ylabel('Density of states (cm$^{-3}$)')
        plt.show()