def display(self, xaxis, alpha, new=True):
        """
        E.display(xaxis, alpha = .8)

        :Arguments: xaxis, alpha

        Plots the CI region on the current figure, with respect to
        xaxis, at opacity alpha.

        :Note: The fill color of the envelope will be self.mass
            on the grayscale.
        """
        if new:
            figure()
        if self.ndim == 1:
            if self.mass>0.:
                x = concatenate((xaxis,xaxis[::-1]))
                y = concatenate((self.lo, self.hi[::-1]))
                fill(x,y,facecolor='%f' % self.mass,alpha=alpha, label = ('centered CI ' + str(self.mass)))
            else:
                pyplot(xaxis,self.value,'k-',alpha=alpha, label = ('median'))
        else:
            if self.mass>0.:
                subplot(1,2,1)
                contourf(xaxis[0],xaxis[1],self.lo,cmap=cm.bone)
                colorbar()
                subplot(1,2,2)
                contourf(xaxis[0],xaxis[1],self.hi,cmap=cm.bone)
                colorbar()
            else:
                contourf(xaxis[0],xaxis[1],self.value,cmap=cm.bone)
                colorbar()

def hinton(W, out_file=None, maxWeight=None):
    """
    Draws a Hinton diagram for visualizing a weight matrix. 
    Temporarily disables matplotlib interactive mode if it is on, 
    otherwise this takes forever.
    """
    reenable = False
    if P.isinteractive():
        P.ioff()
    P.clf()
    height, width = W.shape
    if not maxWeight:
        maxWeight = 2**N.ceil(N.log(N.max(N.abs(W)))/N.log(2))

    P.gca().set_position([0, 0, 1, 1])
    P.fill(N.array([0,width,width,0]),N.array([0,0,height,height]),'gray')
    P.axis('off')
    P.axis('equal')
    for x in xrange(width):
        for y in xrange(height):
            _x = x+1
            _y = y+1
            w = W[y,x]
            if w > 0:
                _blob(_x - 0.5, height - _y + 0.5, min(1,w/maxWeight),'white')
            elif w < 0:
                _blob(_x - 0.5, height - _y + 0.5, min(1,-w/maxWeight),'black')
    if reenable:
        P.ion()
    #P.show()
    if out_file:
        #P.savefig(out_file, format='png', bbox_inches='tight', pad_inches=0)
        fig = P.gcf()
        fig.subplots_adjust()
        fig.savefig(out_file)

def spectrum_subplot (spectrum):
    '''Plot a spectrum, with x-axis the wavelength, and y-axis the intensity.
    The curve is colored at that wavelength by the (approximate) color of a
    pure spectral color at that wavelength, with intensity constant over wavelength.
    (This means that dark looking colors here mean that wavelength is poorly viewed by the eye.

    This is not a complete plotting function, e.g. no file is saved, etc.
    It is assumed that this function is being called by one that handles those things.'''
    (num_wl, num_cols) = spectrum.shape
    # get rgb colors for each wavelength
    rgb_colors = numpy.empty ((num_wl, 3))
    for i in xrange (0, num_wl):
        wl_nm = spectrum [i][0]
        xyz = ciexyz.xyz_from_wavelength (wl_nm)
        rgb_colors [i] = colormodels.rgb_from_xyz (xyz)
    # scale to make brightest rgb value = 1.0
    rgb_max = numpy.max (rgb_colors)
    scaling = 1.0 / rgb_max
    rgb_colors *= scaling        
    # draw color patches (thin vertical lines matching the spectrum curve) in color
    for i in xrange (0, num_wl-1):    # skipping the last one here to stay in range
        x0 = spectrum [i][0]
        x1 = spectrum [i+1][0]
        y0 = spectrum [i][1]
        y1 = spectrum [i+1][1]
        poly_x = [x0,  x1,  x1, x0]
        poly_y = [0.0, 0.0, y1, y0]
        color_string = colormodels.irgb_string_from_rgb (rgb_colors [i])
        pylab.fill (poly_x, poly_y, color_string, edgecolor=color_string)
    # plot intensity as a curve
    pylab.plot (
        spectrum [:,0], spectrum [:,1],
        color='k', linewidth=2.0, antialiased=True)

    def plot_opt_allocation(self, df, file=None):
        try:
            import pylab
        except ImportError:
            pass
        else:
            xs = df[self.symbols].values
            risks = df['std'].values.tolist()
            pylab.figure()

            r = risks[-1::-1] + risks
            cp = [x[0] for x in xs]
            pylab.fill(risks + [.0], cp + [0.0], "g", label=self.symbols[0])

            for i in range(1, len(self.symbols)):
                cn = [sum(x[0:i]) for x in xs]
                c = cn[-1::-1] + cp
                label = self.symbols[i]  # "x%d" % i
                cr = (i * 32 + 0x80) % 255
                cg = (i * 16 + 0x80) % 255
                cb = (i * 64 + 0x80) % 255
                color = "#%02x%02x%02x" % (cr, cg, cb)

                pylab.fill(r, c, label=label, facecolor=color)
                cp = cn

            pylab.legend()
            pylab.xlabel('standard deviation')
            pylab.ylabel('allocation')
            pylab.title('Optimal allocations (fig 4.12)')

            if file is not None:
                pylab.savefig(file)
            else:
                pylab.show()

def plot1d(x,y,cadence,lcolor,lwidth,fcolor,falpha,underfill):

# pad first and last points in case a fill is required

    x = insert(x,[0],[x[0]]) 
    x = append(x,[x[-1]])
    y = insert(y,[0],[-1.0e10]) 
    y = append(y,-1.0e10)

# plot data so that data gaps are not spanned by a line

    ltime = array([],dtype='float64')
    ldata = array([],dtype='float32')
    for i in range(1,len(x)-1):
        if (x[i] - x[i-1]) < 2.0 * cadence / 86400:
            ltime = append(ltime,x[i])
            ldata = append(ldata,y[i])
        else:
            pylab.plot(ltime,ldata,color=lcolor,linestyle='-',linewidth=lwidth)
            ltime = array([],dtype='float64')
            ldata = array([],dtype='float32')
    pylab.plot(ltime,ldata,color=lcolor,linestyle='-',linewidth=lwidth)

# plot the fill color below data time series, with no data gaps

    if underfill:
        pylab.fill(x,y,fc=fcolor,linewidth=0.0,alpha=falpha)

        return

 def fill_gamut_slice (v0, v1, v2):
     '''Fill in a slice of the monitor gamut with the correct colors.'''
     #num_s, num_t = 10, 10
     #num_s, num_t = 25, 25
     num_s, num_t = 50, 50
     dv10 = v1 - v0
     dv21 = v2 - v1
     for i_s in range (num_s):
         s_a = float (i_s)   / float (num_s)
         s_b = float (i_s+1) / float (num_s)
         for i_t in range (num_t):
             t_a = float (i_t)   / float (num_t)
             t_b = float (i_t+1) / float (num_t)
             # vertex coords
             v_aa = v0 + t_a * (dv10 + s_a * dv21)
             v_ab = v0 + t_b * (dv10 + s_a * dv21)
             v_ba = v0 + t_a * (dv10 + s_b * dv21)
             v_bb = v0 + t_b * (dv10 + s_b * dv21)
             # poly coords
             poly_x = [v_aa [0], v_ba [0], v_bb [0], v_ab [0]]
             poly_y = [v_aa [1], v_ba [1], v_bb [1], v_ab [1]]
             # average color
             avg = 0.25 * (v_aa + v_ab + v_ba + v_bb)
             # convert to rgb and scale to maximum displayable brightness
             color_string = get_brightest_irgb_string (avg)
             pylab.fill (poly_x, poly_y, color_string, edgecolor=color_string)

def plot_posterior_ci(locs, mean, sd, color, alpha_multiplier=0.1, rm=True):
    x_ci = SP.array(list(locs) + list(locs)[::-1])
    y_ci = SP.array(list(mean) + list(mean)[::-1])
    if rm: y_ci = 1. - y_ci
    sds = SP.array(list(sd) + list(-sd)[::-1])
    PL.fill(x_ci, y_ci + sds, color, alpha=alpha_multiplier)
    PL.fill(x_ci, y_ci + 2*sds, color, alpha=2*alpha_multiplier) 

def drawDef(dfeat,dy,dx,mindef=0.001,distr="father"):
    """
        auxiliary funtion to draw recursive levels of deformation
    """
    from matplotlib.patches import Ellipse
    pylab.ioff()
    if distr=="father":
        py=[0,0,2,2];px=[0,2,0,2]
    if distr=="child":
        py=[0,1,1,2];px=[1,2,0,1]
    ordy=[0,0,1,1];ordx=[0,1,0,1]
    x1=-0.5+dx;x2=2.5+dx
    y1=-0.5+dy;y2=2.5+dy
    if distr=="father":       
        pylab.fill([x1,x1,x2,x2,x1],[y1,y2,y2,y1,y1],"r", alpha=0.15, edgecolor="b",lw=1)    
    for l in range(len(py)):
        aux=dfeat[ordy[l],ordx[l],:].clip(-1,-mindef)
        wh=numpy.exp(-mindef/aux[0])/numpy.exp(1);hh=numpy.exp(-mindef/aux[1])/numpy.exp(1)
        e=Ellipse(xy=[(px[l]+dx),(py[l]+dy)], width=wh, height=hh, alpha=0.35)
        x1=-0.75+dx+px[l];x2=0.75+dx+px[l]
        y1=-0.76+dy+py[l];y2=0.75+dy+py[l]
        col=numpy.array([wh*hh]*3).clip(0,1)
        if distr=="father":
            col[0]=0       
        e.set_facecolor(col)
        pylab.gca().add_artist(e)
        if distr=="father":       
            pylab.fill([x1,x1,x2,x2,x1],[y1,y2,y2,y1,y1],"b", alpha=0.15, edgecolor="b",lw=1)            

 def myimshow(*args, **kwargs):
     x0,x1,y0,y1 = imExt(afwimg)
     plt.fill([x0,x0,x1,x1,x0],[y0,y1,y1,y0,y0], color=(1,1,0.8),
              zorder=20)
     plt.imshow(*args, zorder=25, **kwargs)
     plt.xticks([]); plt.yticks([])
     plt.axis(imExt(afwimg))

def wiggle(Data,SH,skipt=1,maxval=8,lwidth=.1):
        """
        wiggle(Data,SH)
        """
        import pylab
                
        t = range(SH['ns'])
#       t = range(SH['ns'])*SH['dt']/1000000;

        for i in range(0,SH['ntraces'],skipt):
#               trace=zeros(SH['ns']+2)
#               dtrace=Data[:,i]
#               trace[1:SH['ns']]=Data[:,i]
#               trace[SH['ns']+1]=0
                trace=Data[:,i]
                trace[0]=0
                trace[SH['ns']-1]=0     
                pylab.plot(i+trace/maxval,t,color='black',linewidth=lwidth)
                for a in range(len(trace)):
                        if (trace[a]<0):
                                trace[a]=0;
                # pylab.fill(i+Data[:,i]/maxval,t,color='k',facecolor='g')
                pylab.fill(i+Data[:,i]/maxval,t,'k',linewidth=0)
        pylab.title(SH['filename'])
        pylab.grid(True)

def plot_S2s_over_sequence(S2s_list, label_list, plot_fn, ss_info=None, legend=False, errors_list=None):
    if errors_list != None:
        for S2s, errors, label in zip(S2s_list, errors_list, label_list):
            print range(1, S2s.shape[0]), S2s[1:]
            pylab.errorbar(range(S2s.shape[0]), S2s, fmt="b-", label=label, yerr=errors)
    else:
        for S2s, label in zip(S2s_list, label_list):
            pylab.plot(S2s, "-", label=label)

    # add faded background in SS regions
    if ss_info != None:
        for res_num in ss_info.get_res_nums():
            if ss_info.is_structured(res_num):
                #print "Found SS:", res_num
                x=res_num
                pylab.fill([x-.5,x-.5,x+.5,x+.5], [0,1,1,0], alpha=.3, edgecolor='w')

    #pylab.title("NH order parameters")
    #pylab.ylabel("Order parameter")
    pylab.xlabel("Residue number")
    pylab.ylim(ymax=1)
    pylab.grid()
    if legend: pylab.legend(label_list, prop=matplotlib.font_manager.FontProperties(size='6'), loc='lower right')
    print "Writing ", plot_fn
    pylab.savefig(plot_fn)

def plotSingleYZ(posAll, i):
    ax = pylab.fill([-setup.y-plotFrame,setup.y+plotFrame,setup.y+plotFrame,-setup.y-plotFrame],[-setup.z-plotFrame,-setup.z-plotFrame,setup.z+plotFrame,setup.z+plotFrame],'r')
    bx = pylab.fill([-setup.y,setup.y,setup.y,-setup.y],[-setup.z,-setup.z,setup.z,setup.z],'w')
    pylab.plot([posAll[i][:,1], posAll[i][:,1]], 
               [posAll[i][:,2], posAll[i][:,2]],
               'r.', markersize=5.)
    return

def fibo_boxes(N=10, x_padding=0., y_padding=0., clr_0=None, fill_alpha=.6):
	F=Fibos(N_stop=N)
	plt.figure(0)
	plt.clf()
	#
	colors_ =  mpl.rcParams['axes.color_cycle']
	dy,dx=range(2)
	x=0
	y=0
	for j,f in enumerate(F[1:]):
		side_len=f
		if clr_0==None:
			clr = colors_[j%len(colors_)]
		else:
			clr = clr_0
		#
		square = zip(*[[x,y], [x+side_len, y], [x+side_len,y+side_len], [x, y+side_len], [x,y]])
		print square
		plt.plot(*square, marker='', ls='-', lw=2.5, color=clr)
		plt.fill(*square, color=clr, alpha=fill_alpha)
		#
		x=x+dx*(side_len + x_padding*side_len) - dy*(F[j] + y_padding*side_len)
		y=y+dy*(side_len + y_padding*side_len) - dx*(F[j] + x_padding*side_len)
		
		#
		dx = (1+dx)%2
		dy = (1+dy)%2
	#
	ax=plt.gca()
	ax.set_ylim([-.1*max(square[1]), 1.1*max(square[1])])
	ax.set_xlim([-.1*max(square[0]), 1.1*max(square[0])])

def plotSingleXZ(posAll, i):
    ax = pylab.fill([-setup.x-plotFrame,setup.x+plotFrame,setup.x+plotFrame,-setup.x-plotFrame],[-setup.z-plotFrame,-setup.z-plotFrame,setup.z+plotFrame,setup.z+plotFrame],'r')
    bx = pylab.fill([-setup.x,setup.x,setup.x,-setup.x],[-setup.z,-setup.z,setup.z,setup.z],'w')
    pylab.plot([posAll[i][:,0], posAll[i][:,0]], 
               [posAll[i][:,2], posAll[i][:,2]],
               'k.', markersize=5.)
    return

def plotSingleXY(posAll, i):
    ax = pylab.fill([-setup.x-plotFrame,setup.x+plotFrame,setup.x+plotFrame,-setup.x-plotFrame],[-setup.y-plotFrame,-setup.y-plotFrame,setup.y+plotFrame,setup.y+plotFrame],'r')
    bx = pylab.fill([-setup.x,setup.x,setup.x,-setup.x],[-setup.y,-setup.y,setup.y,setup.y],'w')
    pylab.plot([posAll[i][:,0], posAll[i][:,0]], 
               [posAll[i][:,1], posAll[i][:,1]],
               'k.', markersize=5.)
    return

def hinton(W,filename="hinton.pdf", maxWeight=None):
    """
    Draws a Hinton diagram for visualizing a weight matrix. 
    Temporarily disables matplotlib interactive mode if it is on, 
    otherwise this takes forever.
    """
    reenable = False
    if P.isinteractive():
        P.ioff()
    P.clf()
    ax=P.subplot(111)
    height, width = W.shape
    if not maxWeight:
        maxWeight = 2**np.ceil(np.log(np.max(np.abs(W)))/np.log(2))

    P.fill(np.array([0,width,width,0]),np.array([0,0,height,height]),'gray')
    P.axis('off')
    #P.axis('equal')
    ax.set_yticklabels(['25','20','15','10','5','0'])
    for x in xrange(width):
        for y in xrange(height):
            _x = x+1
            _y = y+1
            w = W[y,x]
            if w > 0:
                _blob(_x - 0.5, height - _y + 0.5, min(1,w/maxWeight),'white')
            elif w < 0:
                _blob(_x - 0.5, height - _y + 0.5, min(1,-w/maxWeight),'black')
    if reenable:
        P.ion()
    P.title(filename)
    P.savefig(filename)

    def plot_spectrum(self, filename, width=0.2, xlim=None):
        """ Make pretty plot of the linear response. 
        
        Parameters:
        ===========
        filename: output file name (&format, supported by matplotlib)
        width:    width of Lorenzian broadening 
        xlim:     energy range for plotting tuple (emin,emax)
        """
        import pylab as pl

        if not self.done:
            self.run()

        e, f = mix.broaden(self.omega * Hartree, self.F, width=width, N=1000, extend=True)
        f = f / max(abs(f))

        pl.plot(e, f, lw=2)
        xs, ys = pl.poly_between(e, 0, f)
        pl.fill(xs, ys, fc="b", ec="b", alpha=0.5)
        pl.ylim(0, 1.2)

        if xlim == None:
            pl.xlim(0, self.emax * Hartree * 1.2)
        else:
            pl.xlim(xlim)
        pl.xlabel("energy (eV)")
        pl.ylabel("linear optical response")
        pl.title("Optical response")
        pl.savefig(filename)
        # pl.show()
        pl.close()

def draw_lorentzian(x,positions,intensities,linewidth,i,vspace):
    for j in range(len(intensities)):
        p = [intensities[j],positions[j],linewidth[j]]
        print p
        lor = Lorentzian(x,p)
        pylab.fill(x,lor+(i*vspace),'b',alpha=0.1)
    return 

    def hinton(W, max_weight=None, names=(names, worst_names)):
        """
        Draws a Hinton diagram for visualizing a weight matrix.
        Temporarily disables matplotlib interactive mode if it is on,
        otherwise this takes forever.
        """
        reenable = False
        if P.isinteractive():
            P.ioff()
        P.clf()
        height, width = W.shape
        if not max_weight:
            max_weight = 2**np.ceil(np.log(np.max(np.abs(W)))/np.log(2))

        P.fill(np.array([0, width, width, 0]), np.array([0, 0, height, height]), 'gray')
        P.axis('off')
        P.axis('equal')
        cmap = plt.get_cmap('RdYlGn')
        for x in range(width):
            if names:
                plt.text(-0.5, x, names[0][x], fontsize=7, ha='right', va='bottom')
                plt.text(x, height+0.5, names[1][height-x-1], fontsize=7, va='bottom', rotation='vertical', ha='left')
            for y in range(height):
                _x = x+1
                _y = y+1
                w = W[y, x]
                if w > 0:
                    _blob(_x - 0.5, height - _y + 0.5, min(1, w/max_weight), color=cmap(w/max_weight))
                elif w < 0:
                    _blob(_x - 0.5, height - _y + 0.5, min(1, -w/max_weight), 'black')
        if reenable:
            P.ion()
        P.show()

    def etch(self, refmask=[0]):
        """Cut away at surface. Can only click polygons for now.
        Optionally input reference mask to guide etching.
        """

        size = self.size
        mask = n.ones((self.size,self.size), dtype='bool')

        print 'Click for points of region to etch (right click to exit).'
        p.figure(1)
        if len(refmask) != 1:
            p.imshow(n.transpose(-refmask), aspect='auto', origin='lower', interpolation='nearest', cmap=p.cm.Greys, extent=(0,1,0,1), vmax=0.5)
        else:
            p.imshow(n.transpose(-self.mask), aspect='auto', origin='lower', interpolation='nearest', cmap=p.cm.Greys, extent=(0,1,0,1), vmax=0.5)
        xy = p.ginput(n=0, timeout=3000)
        xy = n.array(xy)
        
        print 'Calculating remaining region...'
        p.figure(1)
        p.axis([0,1,0,1])
        p.fill(xy[:,0],xy[:,1],'k')

        for i in range(size):
            for j in range(size):
                mask[i,j] = not(point_inside_polygon(float(i)/size, float(j)/size, xy))
        self.mask = self.mask * mask