def plot_preds_orient(weights, plt_num=[1, 1, 1]):
    orients = ['Mean Horiz', 'Mean Vert', 'Mean Diag', 'Mean Diag 2']
    ax = plt.subplot(plt_num[0], plt_num[1], plt_num[2])
    x = 7
    mid = x / 2
    vals_hor = np.zeros([x, x])
    vals_ver = np.zeros([x, x])
    vals_diag = np.zeros([x, x])    
    delta = (mid / 2)
    vals_hor[mid, [mid - delta, mid + delta]] = 1
    vals_ver[[mid - delta, mid + delta], mid] = 1
    x, y = np.diag_indices_from(vals_diag)
    vals_diag[x[mid - delta], y[mid - delta]] = 1
    vals_diag[x[mid + delta], y[mid + delta]] = 1
    vals_diag2 = np.array(zip(*vals_diag[::-1])) 
    
#    delta = mid - mid / 2
#    vals_hor[mid, [mid - delta, mid + delta]] = 1
#    vals_ver[:, mid] = 1
#    np.fill_diagonal(vals_diag, 1)
#    np.fill_diagonal(vals_diag2, 1)
#    vals_diag2 = np.array(zip(*vals_diag2[::-1]))
    
    dat = (vals_hor * weights[0] + vals_ver * weights[1] + vals_diag
           * weights[2] + vals_diag2 * weights[3])
    dat = reverse_fft(dat)
    ext = mid
    _ = plt.imshow(dat, cmap=colors.CoolWarm,
                    interpolation='bilinear',
               extent=(-ext, ext, -ext, ext))
    ylims = np.array([0, np.abs(dat).max()])
    plt.clim(ylims)
    plt.colorbar()
    
    adjust_spines(ax, [])

def run_bootstrap_correlations(Xs, Ys, num_perms, plot=False, verbose=False, permute=False):
# Returns the difference of correlation between X[1]*Y[1] - X[0]*Y[0] using bootstrap
    corrs = []
    pvals = []
    boot_corrs = []

    # need to define ahead of time the bootstrap indexes so we run the same ones for X[1] and X[0] (e.g. baseline and last)
    num_subjects = Xs[0].shape[0]
    boot_idx = np.random.random_integers(0, num_subjects-1, [num_perms, num_subjects])
    for X, Y in zip(Xs, Ys):
        corr = np.empty([X.shape[1], Y.shape[1]])
        pval = np.empty([X.shape[1], Y.shape[1]])
        for x in range(X.shape[1]):
            for y in range(Y.shape[1]):
                corr[x, y], pval[x, y] = stats.pearsonr(X[:, x], Y[:, y])

        corrs.append(corr)
        pvals.append(pval)

        if plot:
            plot_correlations(corr)
            pl.title('n=%d'%X.shape[0])
            pl.clim(-.2, .8)
            pl.draw()

        # checking p-values of the differences
        boot_res = do_bootstrapping(X, Y, num_perms, verbose, reuse_ids=boot_idx)
        boot_corrs.append(boot_res)

    dcorr = boot_corrs[1] - boot_corrs[0]

    return (corrs, pvals, dcorr)

def plot_pairwise_corrcoef(data,ranklist=range(16,24),title="Correlation Coefficient"):
    array_byAttribRankRange=[]
    for attrib in attribList:
        array_byAttribRankRange.append(get_byAttribRankRange(data, attrib=attrib, ranklist=ranklist))
    Narray = len(array_byAttribRankRange)

    array_corrcoef=np.zeros((Narray,Narray),dtype='float')
    for i,elemi in enumerate(array_byAttribRankRange[::-1]):
        for j,elemj in enumerate(array_byAttribRankRange[::-1]):
            if i>j:
                    continue
            elif i==j:
                array_corrcoef[i,j]=1
            else:
                array_corrcoef[i,j]=np.corrcoef(elemi,elemj)[0,1]

    P.pcolor(np.transpose(array_corrcoef), cmap=P.cm.RdBu, alpha=0.8)
    P.title(title)
    P.xlim([0,23])
    P.ylim([0,23])
    P.clim([-1,1])
    P.xticks(range(len(attribList)), attribListAbrv[::-1],rotation='vertical')
    P.yticks(range(len(attribList)), attribListAbrv[::-1])
    P.subplots_adjust(bottom=0.35)
    P.subplots_adjust(left=0.25)
    P.colorbar()
    return array_corrcoef

    def __call__(self, n):
        if len(self.f.shape) == 3:
            # f = f[x,v,t], 2 dim in phase space
            ft = self.f[n,:,:]
            pylab.pcolormesh(self.X, self.V, ft.T, cmap = 'jet')
            pylab.colorbar()
            pylab.clim(0,0.38) # for Landau test case
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$v$', fontsize = 18)
            pylab.title('$N_x$ = %d, $N_v$ = %d, $t$ = %2.1f' % (self.x.N, self.v.N, self.it*self.t.width))
            pylab.savefig(self.path + self.filename)
            pylab.clf()
            return None

        if len(self.f.shape) == 2:
            # f = f[x], 1 dim in phase space
            ft = self.f[n,:]
            pylab.plot(self.x.gridvalues,ft,'ob')
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$f(x)$', fontsize = 18)
            pylab.savefig(self.path + self.filename)
            return None

    def plot_onereg(self, regid):
        """Plot mintime of connectivity for one individual region"""
        tomat = self.llat * np.nan
        for i,v in enumerate(self.mintmat[1:,regid]):
            tomat[self.regmat==(i+1)] = v
        frmat = self.llat * np.nan
        for i,v in enumerate(self.mintmat[regid,1:]):
            frmat[self.regmat==(i+1)] = v
        djtk,lim = djdticks(max(np.nanmax(tomat), np.nanmax(frmat)))
        figpref.current()
        pl.close(1)
        fig = pl.figure(1,(8,8))
        pl.clf()
        pl.suptitle("Connectivity for region %i" % regid)

        mask = self.regmat == regid
        x,y = self.gcm.mp(self.llon[mask], self.llat[mask])
        pl.subplots_adjust(hspace=0, top=0.95,bottom = 0.15)

        pl.subplot(2,1,1, axisbg="0.8")
        self.gcm.pcolor(frmat, cmap=WRY(), rasterized=True)
        pl.clim(0,lim)
        self.gcm.mp.text(70,60,'Time from region')
        self.gcm.mp.scatter(x,y,5,'b')

        pl.subplot(2,1,2, axisbg="0.8")
        self.gcm.pcolor(tomat, cmap=WRY(), rasterized=True)
        pl.clim(0,lim)
        self.gcm.mp.text(70,60,'Time to region')
        self.gcm.mp.scatter(x,y,5,'b')

        if 'mycolor' in sys.modules:
            mycolor.freecbar([0.2,0.12,0.6,0.025], djtk, cmap=WRY())
        pl.savefig('figs/onereg_%02i_%02i_%06i.png' %
                   (self.regdi, self.regdj, regid))

def transect(x,y,z,x0,y0,x1,y1,plots=0):
    #convert coord to pixel coord
    d0=sqrt( (x-x0)**2+ (y-y0)**2 );
    i0=d0.argmin();
    x0,y0=unravel_index(i0,x.shape); #overwrite x0,y0
    
    d1=plt.np.sqrt( (x-x1)**2+ (y-y1)**2 );
    i1=d1.argmin();
    x1,y1=unravel_index(i1,x.shape); #overwrite x1,y1    
    #-- Extract the line...    
    # Make a line with "num" points...
    length = int(plt.np.hypot(x1-x0, y1-y0))
    xi, yi = plt.np.linspace(x0, x1, length), plt.np.linspace(y0, y1, length) 
       
    # Extract the values along the line
    #y is the first dimension and x is the second, row,col
    zi = z[xi.astype(plt.np.int), yi.astype(plt.np.int)]
    mz=nonaninf(z.ravel()).mean()
    sz=nonaninf(z.ravel()).std()
    if plots==1:
        plt.matshow(z);plt.clim([mz-2*sz,mz+2*sz]);plt.colorbar();plt.title('transect: (' + str(x0) + ',' + str(y0) + ') (' +str(x1) + ',' +str(y1) + ')' );
        plt.scatter(yi,xi,5,c='r',edgecolors='none')
        plt.figure();plt.scatter(sqrt( (xi-xi[0])**2 + (yi-yi[0])**2 ) , zi)
        #plt.figure();plt.scatter(xi, zi)
        #plt.figure();plt.scatter(yi, zi)

    return (xi, yi, zi);

def show_kappa(kappa,clabel=r'\kappa',part='r',clim=None,logplot=True):
    if part.lower()=='r':
        kappa = kappa.real
    elif part.lower()=='i':
        kappa = kappa.imag
    elif part.lower() in ['a','n']:
        kappa = abs(kappa)
    else:
        raise ValueError, "show_kappa : unrecognized part %s" % part
    
    pylab.figure(figsize=(14,9))

    if logplot:
        kappa = numpy.log(1+kappa)
    
    pylab.imshow(kappa.T,
                 origin = 'lower',
                 interpolation = 'nearest',
                 extent = params.RAlim+params.DEClim)
    cb = pylab.colorbar()
    if logplot:
        cb.set_label(r'$\rm{log}(1+%s)$' % clabel,
                     fontsize=14)
    else:
        cb.set_label(r'$%s$' % clabel,
                     fontsize=14)

    if clim is not None:
        pylab.clim(clim)
    
    pylab.xlabel('RA (deg)')
    pylab.ylabel('DEC (deg)')

 def movie(self):
     import matplotlib as mpl
     mpl.rcParams['axes.labelcolor'] = 'white'
     pl.close(1)
     pl.figure(1,(8,4.5),facecolor='k')
     miv = np.ma.masked_invalid
     figpref.current()
     jd0 = pl.date2num(dtm(2005,1,1))
     jd1 = pl.date2num(dtm(2005,12,31))
     mp = projmaps.Projmap('glob')
     x,y = mp(self.llon,self.llat)
     for t in np.arange(jd0,jd1):
         print pl.num2date(t)
         self.load(t)
     
         pl.clf()
         pl.subplot(111,axisbg='k')
         mp.pcolormesh(x,y,
                       miv(np.sqrt(self.u**2 +self.v**2)),
                       cmap=cm.gist_heat)
         pl.clim(0,1.5)
         mp.nice()
         pl.title('%04i-%02i-%02i' % (pl.num2date(t).year,
                                      pl.num2date(t).month,
                                      pl.num2date(t).day),
                  color='w')
         pl.savefig('/Users/bror/oscar/norm/%03i.png' % t,
                    bbox_inches='tight',facecolor='k',dpi=150)

def showqtresultfit(thk, wc, t2, datvec, resp, t,
                    islog=True, clim=None, nu=3, nv=2):
    ''' show mrs qt result and data fit
    showqtresultfit(thk,wc,t2,datvec,resp,t,islog=True,clim=None,nu=3,nv=2)
    '''
    if clim is None:
        cma = max(datvec)
        cmi = min(datvec)
        if islog:
            cma = N.log10(cma)
            cmi = cma - 1.5
        clim = (cmi, cma)

    nt = len(t)
    nq = len(datvec) / nt
    si = (nq, nt)

#    P.clf()
#    P.subplot(nu,nv,1)

    fig = P.figure(1)
    ax1 = fig.add_subplot(nu, nv, 1)

    draw1dmodel(wc, thk, islog=False, xlab=r'$\theta$')
#    P.subplot(nu,nv,3)
    ax3 = fig.add_subplot(nu, nv, 3)
    draw1dmodel(t2, thk, xlab='T2* in ms')
    ax3.set_xticks( [0.02, 0.05, 0.1, 0.2, 0.5] )
    ax3.set_xticklabels( ('0.02', '0.05', '0.1', '0.2', '0.5') )
#    P.subplot(nu,nv,2)
    ax2 = fig.add_subplot(nu, nv, 2)
    if islog:
        P.imshow(N.log10( N.array(datvec).reshape(si) ),
                 interpolation='nearest', aspect='auto')
    else:
        P.imshow(N.array(datvec).reshape(si),
                 interpolation='nearest', aspect='auto')
    P.clim(clim)
#    P.subplot(nu,nv,4)
    ax4 = fig.add_subplot(nu, nv, 4)
    if islog:
        P.imshow(N.log10(resp.reshape(si)),
                 interpolation='nearest',aspect='auto')
    else:
        P.imshow(resp.reshape(si),
                 interpolation='nearest',aspect='auto')
    misfit = N.array( datvec - resp )
    P.clim(clim)
#    P.subplot(nu,nv,5)
    ax5 = fig.add_subplot(nu, nv, 5)
    P.hist(misfit, bins=30)
    P.axis('tight')
    P.grid(which='both')
    P.text(P.xlim()[0], N.mean( P.ylim() ),
           ' std=%g nV' % rndig( N.std(misfit), 3 ) )
#    P.subplot(nu,nv,6)
    ax6 = fig.add_subplot(nu, nv, 6)
    P.imshow(misfit.reshape(si), interpolation='nearest', aspect='auto')
    ax = [ ax1, ax2, ax3, ax4, ax5, ax6 ]
    return ax

def three_component_plot(c1, c2, c3, color, labels):
    pl.figure(figsize=(8,8))
    kwargs = dict(s=4, lw=0, c=color, vmin=2, vmax=6)
    ax1 = pl.subplot(221)
    pl.scatter(c1, c2, **kwargs)
    pl.ylabel('component 2')

    ax2 = pl.subplot(223, sharex=ax1)
    pl.scatter(c1, c3, **kwargs)
    pl.xlabel('component 1')
    pl.ylabel('component 3')

    ax3 = pl.subplot(224, sharey=ax2)
    pl.scatter(c2, c3, **kwargs)
    pl.xlabel('component 2')

    for ax in (ax1, ax2, ax3):
        ax.xaxis.set_major_formatter(ticker.NullFormatter())
        ax.yaxis.set_major_formatter(ticker.NullFormatter())

    pl.subplots_adjust(hspace=0.05, wspace=0.05)

    format = ticker.FuncFormatter(lambda i, *args: labels[i])
    pl.colorbar(ticks = range(2, 7), format=format,
                cax = pl.axes((0.52, 0.51, 0.02, 0.39)))
    pl.clim(1.5, 6.5)

def plot_modes_6(filename,modes=[0,1,2,3,49,499]):
    assert(len(modes)==6) 
    X = numpy.load(filename)
    evecs = X['evecs']
    print evecs.shape

    N = int( numpy.sqrt(evecs.shape[0]) )

    pylab.figure(figsize=(7,10))
    for i in range(6):
        evec = evecs[:,modes[i]]
        pylab.subplot(321+i)
        pylab.imshow(evec.reshape((N,N)),
                     origin='lower',
                     cmap=pylab.cm.RdGy,
                     extent=(0,60,0,60) )

        pylab.title('n=%i' % (modes[i]+1))
        pylab.colorbar()
        cmax = numpy.max(abs(evec))
        pylab.clim(-cmax,cmax)
        pylab.xlabel('arcmin')
        
        if i%2 == 1:
            pylab.gca().yaxis.set_major_formatter(NullFormatter())
        else:
            pylab.ylabel('arcmin')

def plot_color_map_2d( data,  interpolation = None, limits = None, 
                       axis_labels = None, title = None, show = False ):
    """
    interpolation = None will default to using interpolation
    interpolation = 'none' will ensure no interpolation is performed
    """
    assert data.ndim == 2
    p = pylab.imshow( data, interpolation = interpolation )
    fig = pylab.gcf()
    pylab.clim()   # clamp the color limits
    pylab.colorbar() # add color bar to indicate scale
    
    if limits is not None:
        xlim = limits[:2]
        ylim = limits[2:]

        num_ticks = 11

        M, N = data.shape
        xticks = numpy.linspace( 0, M-1, num_ticks )
        yticks = numpy.linspace( 0, N-1, num_ticks )
        
        pylab.xticks( xticks, numpy.linspace( xlim[0], xlim[1], num_ticks ) )
        pylab.yticks( yticks, numpy.linspace( ylim[0], ylim[1], num_ticks ) )

    if axis_labels is not None:
        p.set_xlabel( axis_labels[0] )
        p.set_ylabel( axis_labels[1] )

    if title is not None:
        pylab.title( title )

    if show:
        pylab.show()

def plotWeightChanges():
    if f.usestdp:
        # create plot
        figh = figure(figsize=(1.2*8,1.2*6))
        figh.subplots_adjust(left=0.02) # Less space on left
        figh.subplots_adjust(right=0.98) # Less space on right
        figh.subplots_adjust(top=0.96) # Less space on bottom
        figh.subplots_adjust(bottom=0.02) # Less space on bottom
        figh.subplots_adjust(wspace=0) # More space between
        figh.subplots_adjust(hspace=0) # More space between
        h = axes()

        # create data matrix
        wcs = [x[-1][-1] for x in f.allweightchanges] # absolute final weight
        wcs = [x[-1][-1]-x[0][-1] for x in f.allweightchanges] # absolute weight change
        pre,post,recep = zip(*[(x[0],x[1],x[2]) for x in f.allstdpconndata])
        ncells = int(max(max(pre),max(post))+1)
        wcmat = zeros([ncells, ncells])

        for iwc,ipre,ipost,irecep in zip(wcs,pre,post,recep):
            wcmat[int(ipre),int(ipost)] = iwc *(-1 if irecep>=2 else 1)

        # plot
        imshow(wcmat,interpolation='nearest',cmap=bicolormap(gap=0,mingreen=0.2,redbluemix=0.1,epsilon=0.01))
        xlabel('post-synaptic cell id')
        ylabel('pre-synaptic cell id')
        h.set_xticks(f.popGidStart)
        h.set_yticks(f.popGidStart)
        h.set_xticklabels(f.popnames)
        h.set_yticklabels(f.popnames)
        h.xaxif.set_ticks_position('top')
        xlim(-0.5,ncells-0.5)
        ylim(ncells-0.5,-0.5)
        clim(-abs(wcmat).max(),abs(wcmat).max())
        colorbar()

def plotConn():
    # Create plot
    figh = figure(figsize=(8,6))
    figh.subplots_adjust(left=0.02) # Less space on left
    figh.subplots_adjust(right=0.98) # Less space on right
    figh.subplots_adjust(top=0.96) # Less space on bottom
    figh.subplots_adjust(bottom=0.02) # Less space on bottom
    figh.subplots_adjust(wspace=0) # More space between
    figh.subplots_adjust(hspace=0) # More space between
    h = axes()
    totalconns = zeros(shape(f.connprobs))
    for c1 in range(size(f.connprobs,0)):
        for c2 in range(size(f.connprobs,1)):
            for w in range(f.nreceptors):
                totalconns[c1,c2] += f.connprobs[c1,c2]*f.connweights[c1,c2,w]*(-1 if w>=2 else 1)
    imshow(totalconns,interpolation='nearest',cmap=bicolormap(gap=0))

    # Plot grid lines
    hold(True)
    for pop in range(f.npops):
        plot(array([0,f.npops])-0.5,array([pop,pop])-0.5,'-',c=(0.7,0.7,0.7))
        plot(array([pop,pop])-0.5,array([0,f.npops])-0.5,'-',c=(0.7,0.7,0.7))

    # Make pretty
    h.set_xticks(range(f.npops))
    h.set_yticks(range(f.npops))
    h.set_xticklabels(f.popnames)
    h.set_yticklabels(f.popnames)
    h.xaxis.set_ticks_position('top')
    xlim(-0.5,f.npops-0.5)
    ylim(f.npops-0.5,-0.5)
    clim(-abs(totalconns).max(),abs(totalconns).max())
    colorbar()

def plot_eigenmodes(evecs, evals, ax_list, mode_list, RArange, DECrange):
    """
    Plot KL eigenmodes associated with the COSMOS catalog
    """
    assert len(ax_list) == len(mode_list)

    NRA = len(RArange) - 1
    NDEC = len(DECrange) - 1

    for i in range(len(ax_list)):
        ax = ax_list[i]
        mode = mode_list[i]

        pylab.axes(ax)

        evec = evecs[:, i]

        pylab.imshow(
            evec.reshape((NRA, NDEC)).T,
            origin="lower",
            interpolation=None,  #'nearest',
            cmap=pylab.cm.RdGy,
            extent=(RArange[0], RArange[-1], DECrange[0], DECrange[-1]),
        )
        cmax = np.max(abs(evec))
        pylab.clim(-cmax, cmax)

        pylab.title(r"$\mathrm{mode\ %i}\ (v=%.3f)$" % (i + 1, evals[i]))

    return ax_list

 def plot_discharge(self,grid):
     """This plots discharge across the raster"""
  
     grid._setup_active_inlink_and_outlink_matrices()
     outlink = grid.node_active_outlink_matrix
     inlink = grid.node_active_inlink_matrix
     outlink1 = outlink.tolist()
     inlink1 =inlink.tolist()
     newin0, newout0 = change_signs(inlink1[0], outlink1[0])
     newin1, newout1 = change_signs(inlink1[1], outlink1[1])
     in0 = np.array(newin0)
     in1 = np.array(newin1)
     out0 = np.array(newout0)
     out1 = np.array(newout1)
     self.q_node = self.q[in0]+self.q[in1]+self.q[out0]+self.q[out1] #((q[outlink[1]] -q[inlink[1]]))#+((q[outlink[0]]-(q[inlink[0]])))
     fixed_q = grid.zeros(centering='node')
     for each in self.interior_nodes:
         fixed_q[each] = self.q_node[each]
     plt.figure('DISCHARGE')
     hr = grid.node_vector_to_raster(fixed_q)
     palette = pylab.cm.RdYlBu
     im2 = pylab.imshow(hr, cmap=palette, extent=[0, grid.number_of_node_columns *grid.dx,0, grid.number_of_node_rows * grid.dx])
     pylab.clim(vmin=0.000001)#, mx)
     palette.set_under('w', 0.000001)
     cb = pylab.colorbar(im2)
     cb.set_label('DISCHARGE (m)', fontsize=12)
     pylab.title('DISCHARGE')
     plt.show()

def plot_mtx(mtx, labels, title):
    pl.figure()
    pl.imshow(mtx, interpolation='nearest')
    pl.xticks(range(len(mtx)), labels, rotation=-45)
    pl.yticks(range(len(mtx)), labels)
    pl.title(title)
    pl.clim((0,1))
    pl.colorbar()

def saltmovie(tr):
    mv = anim.Movie()
    ints = np.sort(np.unique(tr.ints))
    for i in ints[:40]:
        pl.pcolormesh(miv(tr.traj2grid('z',tr.jd==i)),cmap=cm.Paired)
        pl.clim(28,33)
        print i
        mv.image()
    mv.video(tr.projname+tr.casename+"_salt_mov.mp4")

 def animate(i):
     im.set_array(data[i])
     if saturation is not None:
         saturate(saturation, im=im)
     else:
         plt.clim(clims)
     if label is not None:
         tx.set_text('{:2.1f} {:}'.format(label[i], label_caption))
     return [im]

    def movie(self,fld,k=0,jd1=733342,jd2=733342+10):
        miv = np.ma.masked_invalid

        for n,jd in enumerate(np.arange(jd1,jd2,0.25)):
            self.load(fld,jd)
            pl.clf()
            pl.pcolormesh(miv(self.__dict__[fld][0,k,:,:]))
            pl.clim(-0.1,0.1)
            pl.savefig('%s_%05i.png' % (fld,n),dpi=150)