def window_fn_matrix(Q,N,num_remov=None,save_tag=None,lms=None):
    Q = n.matrix(Q); N = n.matrix(N)
    Ninv = uf.pseudo_inverse(N,num_remov=None) # XXX want to remove dynamically
    #print Ninv 
    info = n.dot(Q.H,n.dot(Ninv,Q))
    M = uf.pseudo_inverse(info,num_remov=num_remov)
    W = n.dot(M,info)

    if save_tag!=None:
        foo = W[0,:]
        foo = n.real(n.array(foo))
        foo.shape = (foo.shape[1]),
        print foo.shape
        p.scatter(lms[:,0],foo,c=lms[:,1],cmap=mpl.cm.PiYG,s=50)
        p.xlabel('l (color is m)')
        p.ylabel('W_0,lm')
        p.title('First Row of Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W.pdf'.format(fig_loc,save_tag))
        p.clf()

        print 'W ',W.shape
        p.imshow(n.real(W))
        p.title('Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W_im.pdf'.format(fig_loc,save_tag))
        p.clf()


    return W

def plotear(xi,yi,zi):
    # mask inner circle
    interior1 = sqrt(((xi+1.5)**2) + (yi**2)) < 1.0 
    interior2 = sqrt(((xi-1.5)**2) + (yi**2)) < 1.0
    zi[interior1] = ma.masked
    zi[interior2] = ma.masked
    p.figure(figsize=(16,10))
    pyplot.jet()
    max=2.8
    min=0.4
    steps = 50
    levels=list()
    labels=list()
    for i in range(0,steps):
	levels.append(int((max-min)/steps*100*i)*0.01+min)
    for i in range(0,steps/2):
	labels.append(levels[2*i])
    CSF = p.contourf(xi,yi,zi,levels,norm=colors.LogNorm())
    CS = p.contour(xi,yi,zi,levels, format='%.3f', labelsize='18')
    p.clabel(CS,labels,inline=1,fontsize=9)
    p.title('electrostatic potential of two spherical colloids, R=lambda/3',fontsize=24)
    p.xlabel('z-coordinate (3*lambda)',fontsize=18)
    p.ylabel('radial coordinate r (3*lambda)',fontsize=18)
    # add a vertical bar with the color values
    cbar = p.colorbar(CSF,ticks=labels,format='%.3f')
    cbar.ax.set_ylabel('potential (reduced units)',fontsize=18)
    cbar.add_lines(CS)
    p.show()

    def Xtest3(self):
        """
        Test from Kate Marvel
        As the following code snippet demonstrates, regridding a
        cdms2.tvariable.TransientVariable instance using regridTool='regrid2' 
        results in a new array that is masked everywhere.  regridTool='esmf' 
        and regridTool='libcf' both work as expected.

        This is similar to the original test but we construct our own 
        uniform grid. This should passes.
        """
        import cdms2 as cdms
        import numpy as np

        filename = cdat_info.get_sampledata_path() + '/clt.nc'
        a=cdms.open(filename)
        data=a('clt')[0,...]

        print data.mask #verify this data is not masked

        GRID = cdms.grid.createUniformGrid(-90.0, 23, 8.0, -180.0, 36, 10.0, order="yx", mask=None)

        test_data=data.regrid(GRID,regridTool='regrid2')

        # check that the mask does not extend everywhere...
        self.assertNotEqual(test_data.mask.sum(), test_data.size)
        
        if PLOT:
            pylab.subplot(2, 1, 1)
            pylab.pcolor(data[...])
            pylab.title('data')
            pylab.subplot(2, 1, 2)
            pylab.pcolor(test_data[...])
            pylab.title('test_data (interpolated data)')
            pylab.show()

	def plot_stress(self, block_ids=None, fignum=0):
		block_ids = self.check_block_ids_list(block_ids)
		#
		plt.figure(fignum)
		ax1=plt.gca()
		plt.clf()
		plt.figure(fignum)
		plt.clf()
		ax0=plt.gca()

		#
		for block_id in block_ids:
			rws = numpy.core.records.fromarrays(zip(*filter(lambda x: x['block_id']==block_id, self.shear_stress_sequences)), dtype=self.shear_stress_sequences.dtype)
			stress_seq = []
			for rw in rws:
				stress_seq += [[rw['sweep_number'], rw['shear_init']]]
				stress_seq += [[rw['sweep_number'], rw['shear_final']]]
			X,Y = zip(*stress_seq)
			#
			ax0.plot(X,Y, '.-', label='block_id: %d' % block_id)
			#
			plt.figure(fignum+1)
			plt.plot(rws['sweep_number'], rws['shear_init'], '.-', label='block_id: %d' % block_id)
			plt.plot(rws['sweep_number'], rws['shear_final'], '.-', label='block_id: %d' % block_id)
			plt.figure(fignum)
		ax0.plot([min(self.shear_stress_sequences['sweep_number']), max(self.shear_stress_sequences['sweep_number'])], [0., 0.], 'k-')
		ax0.legend(loc=0, numpoints=1)
		plt.figure(fignum)
		plt.title('Block shear_stress sequences')
		plt.xlabel('sweep number')
		plt.ylabel('shear stress')

 def showHistory(self, figNum):
     pylab.figure(figNum)
     plot = pylab.plot(self.history, label = 'Test Stock')
     plot
     pylab.title('Closing Price, Test ' + str(figNum))
     pylab.xlabel('Day')
     pylab.ylabel('Price')

    def add_quad(self, name, X, fZ):
        """
        Create a 4-panel figure
        name: title of the figure
        X: 1D axes data
        fZ: 2D data set (Fortran indexing)
        """
# Does not work on chipolata: (older matplotlib?)
#	i = 0
#	fig, axs = pylab.subplots(2,2)
#	for ax in axs.ravel():
#		pcm = ax.pcolormesh(F[i])
#		ax.axis('tight')
#		fig.colorbar(pcm,ax=ax)
#		ax.set_title(name+','+str(i))
#		i += 1

# The default ordering for 2D meshgrid is Fortran style
        title = name.replace(',', '_')
        self.figs.append((title, pylab.figure(title)))
        for i in range(4):
            fX, fY = numpy.meshgrid(X[i], X[i])
            ax = pylab.subplot('22'+str(i+1))
            pcm = ax.pcolormesh(fX, fY, fZ[i])
            ax.axis('tight')
            self.figs[-1][1].colorbar(pcm)
            pylab.title(name+','+str(i))

def plotHistogram(data, preTime):
    pylab.figure(1)
    pylab.hist(data, bins=10)
    pylab.xlabel("Virus Population At End of Simulation")
    pylab.ylabel("Number of Trials")
    pylab.title("{0} Time Steps Before Treatment Simulation".format(preTime))
    pylab.show()

def plotB3reg():
    w=loadStanFit('revE2B3BHreg.fit')
    printCI(w,'mmu')
    printCI(w,'mr')
    for b in range(2):
        subplot(1,2,b+1)
        plt.title('')
        px=np.array(np.linspace(-0.5,0.5,101),ndmin=2)
        a0=np.array(w['mmu'][:,b],ndmin=2).T
        a1=np.array(w['mr'][:,b],ndmin=2).T
        y=np.concatenate([sap(a0+a1*px,97.5,axis=0),sap(a0+a1*px[:,::-1],2.5,axis=0)])
        x=np.squeeze(np.concatenate([px,px[:,::-1]],axis=1))
        plt.plot(px[0,:],np.median(a0)+np.median(a1)*px[0,:],'red')
        #plt.plot([-1,1],[0.5,0.5],'grey')
        ax=plt.gca()
        ax.set_aspect(1)
        ax.add_patch(plt.Polygon(np.array([x,y]).T,alpha=0.2,fill=True,fc='red',ec='w'))
        y=np.concatenate([sap(a0+a1*px,75,axis=0),sap(a0+a1*px[:,::-1],25,axis=0)])
        ax.add_patch(plt.Polygon(np.array([x,y]).T,alpha=0.2,fill=True,fc='red',ec='w'))
        man=np.array([-0.4,-0.2,0,0.2,0.4])
        mus=[]
        for m in range(len(man)):
            mus.append(loadStanFit('revE2B3BH%d.fit'%m)['mmu'][:,b])
        mus=np.array(mus).T
        errorbar(mus,x=man)
        ax.set_xticks(man)
        plt.xlim([-0.5,0.5])
        plt.ylim([-0.4,0.8])
        #plt.xlabel('Manipulated Displacement')
        if b==0:
            plt.ylabel('Perceived Displacemet')
            plt.gca().set_yticklabels([])
        subplot_annotate()
    plt.text(-1.1,-0.6,'Pivot Displacement',fontsize=8);

    def plot_cost(self):
        if self.show_cost not in self.train_outputs[0][0]:
            raise ShowNetError("Cost function with name '%s' not defined by given convnet." % self.show_cost)
        train_errors = [o[0][self.show_cost][self.cost_idx] for o in self.train_outputs]
        test_errors = [o[0][self.show_cost][self.cost_idx] for o in self.test_outputs]

        numbatches = len(self.train_batch_range)
        test_errors = numpy.row_stack(test_errors)
        test_errors = numpy.tile(test_errors, (1, self.testing_freq))
        test_errors = list(test_errors.flatten())
        test_errors += [test_errors[-1]] * max(0,len(train_errors) - len(test_errors))
        test_errors = test_errors[:len(train_errors)]

        numepochs = len(train_errors) / float(numbatches)
        pl.figure(1)
        x = range(0, len(train_errors))
        pl.plot(x, train_errors, 'k-', label='Training set')
        pl.plot(x, test_errors, 'r-', label='Test set')
        pl.legend()
        ticklocs = range(numbatches, len(train_errors) - len(train_errors) % numbatches + 1, numbatches)
        epoch_label_gran = int(ceil(numepochs / 20.)) # aim for about 20 labels
        epoch_label_gran = int(ceil(float(epoch_label_gran) / 10) * 10) # but round to nearest 10
        ticklabels = map(lambda x: str((x[1] / numbatches)) if x[0] % epoch_label_gran == epoch_label_gran-1 else '', enumerate(ticklocs))

        pl.xticks(ticklocs, ticklabels)
        pl.xlabel('Epoch')
#        pl.ylabel(self.show_cost)
        pl.title(self.show_cost)

def evaluate_result(data, target, result):
	assert(data.shape[0] == target.shape[0])
	assert(target.shape[0] == result.shape[0])
	
	correct = np.where( result == target )
	miss 	= np.where( result != target )
	
	class_rate = float(correct[0].shape[0]) / target.shape[0]

	print "Correct classification rate:", class_rate 
	#get the 3s
	mask 			= np.where(target == wanted[0])
	data_3_correct 	= data[np.intersect1d(mask[0],correct[0])]
	data_3_miss	 	= data[np.intersect1d(mask[0],miss[0])]
	#get the 8s
	mask = np.where(target == wanted[1])
	data_8_correct 	= data[np.intersect1d(mask[0],correct[0])]
	data_8_miss	 	= data[np.intersect1d(mask[0],miss[0])]
	#plot
	plot.title("Scatter")
	plot.xlabel("x_0")
	plot.ylabel("x_1")
	size = 20
	plot.scatter(data_3_correct[:,0], data_3_correct[:,1], marker = "x", c = "r", s = size )
	plot.scatter(   data_3_miss[:,0],    data_3_miss[:,1], marker = "x", c = "b", s = size )
	plot.scatter(data_8_correct[:,0], data_8_correct[:,1], marker = "o", c = "r", s = size )
	plot.scatter(   data_8_miss[:,0],    data_8_miss[:,1], marker = "o", c = "b", s = size )
	plot.show()

    def plot_roc(self, roc=None):
        """Plot ROC curves

        .. plot::
            :include-source:
            :width: 80%

            from dreamtools import rocs
            r = rocs.ROC()
            r.scores = [.9,.5,.6,.7,.1,.2,.6,.4,.7,.9, .2]
            r.classes = [1,0,1,0,0,1,1,0,0,1,1]
            r.plot_roc()

        """
        if roc == None:
            roc = self.get_roc()
        from pylab import plot, xlim, ylim ,grid, title, xlabel, ylabel
        x = roc['fpr']
        plot(x, roc['tpr'], '-o')
        plot([0,1], [0,1],'r')
        ylim([0, 1])
        xlim([0, 1])
        grid(True)
        title("ROC curve (AUC=%s)" % self.compute_auc(roc))
        xlabel("FPR")
        ylabel("TPR")

def PlotLine(type):
	i=j=0
	pylab.figure(type)
	if type==0:
		pylab.title("Geomagnetism")
		pylab.xlabel("Distance")
		pylab.ylabel("Value")
	elif type==1:
		pylab.title("Compass")
		pylab.xlabel("Distance")
		pylab.ylabel("Value")
	for path in pathlist:
		f = open(path)
		f.readline()
		data = np.loadtxt(f)
		dataAfterfilter = filters.median(data,10)
		if type == 0:
			pylab.plot(dataAfterfilter,  color[i] ,label =lablelist[i])
			i=i+1
			pylab.legend()
		elif type == 1:
			pylab.plot(data[:,1],  color[i] ,label =lablelist[i])
			i=i+1
			pylab.legend()
	pass

def simulationWithDrug():

    """

    Runs simulations and plots graphs for problem 4.
    Instantiates a patient, runs a simulation for 150 timesteps, adds
    guttagonol, and runs the simulation for an additional 150 timesteps.
    total virus population vs. time and guttagonol-resistant virus population
    vs. time are plotted
    """

    maxBirthProb = .1
    clearProb = .05
    resistances = {'guttagonal': False}
    mutProb = .005

    total = [100]
    g = [0]

    badVirus = ResistantVirus(maxBirthProb, clearProb, resistances, mutProb)
    viruses = [badVirus]*total[0]
    maxPop = 1000

    Bob = Patient(viruses, maxPop)

    for i in range(150):
        Bob.update()
        gVirus = 0

        for v in Bob.viruses:
            if v.isResistantTo('guttagonal'):
                gVirus += 1

        #print "g = ", gVirus
        #print "t = ", len(Bob.viruses)
        #print
        g += [gVirus]
        total += [len(Bob.viruses)]

    Bob.addPrescription('guttagonal')

    for i in range(150):
        Bob.update()
        gVirus = 0

        for v in Bob.viruses:
            if v.isResistantTo('guttagonal'):
                gVirus += 1

        g += [gVirus]
        total += [len(Bob.viruses)]

    pylab.title("Number of Viruses with Different Resistances to Guttagonal")
    pylab.xlabel("Number of Timesteps")
    pylab.ylabel("Number of Viruses")

    pylab.plot(g, '-r', label = 'Resistant')
    pylab.plot(total, '-b', label = 'Total')
    pylab.legend(loc = 'lower right')
    pylab.show()

def plot_signal(x,y,title,labelx,labley,position):
    pylab.subplot (9, 1, position)
    pylab.plot(x,y)
    pylab.title(title)
    pylab.xlabel(labelx)
    pylab.ylabel(labley)
    pylab.grid(True)

def drawPrfastscore(tp,fp,scr,tot,show=True):
    tp=numpy.cumsum(tp)
    fp=numpy.cumsum(fp)
    rec=tp/tot
    prec=tp/(fp+tp)
    #dif=numpy.abs(prec[1:]-rec[1:])
    dif=numpy.abs(prec[::-1]-rec[::-1])
    pos=dif.argmin()
    pos=len(dif)-pos-1
    ap=0
    for t in numpy.linspace(0,1,11):
        pr=prec[rec>=t]
        if pr.size==0:
            pr=0
        p=numpy.max(pr);
        ap=ap+p/11;
    if show:    
        pylab.plot(rec,prec,'-g')
        pylab.title("AP=%.3f EPRthr=%.3f"%(ap,scr[pos]))
        pylab.xlabel("Recall")
        pylab.ylabel("Precision")
        pylab.grid()
        pylab.show()
        pylab.draw()
    return rec,prec,scr,ap,scr[pos]

def Doplots_monthly(mypathforResults,PlottingDF,variable_to_fill, Site_ID,units,item):   
    ANN_label=str(item+"_NN")     #Do Monthly Plots
    print "Doing MOnthly  plot"
    #t = arange(1, 54, 1)
    NN_label='Fc'
    Plottemp = PlottingDF[[NN_label,item]][PlottingDF['day_night']!=1]
    #Plottemp = PlottingDF[[NN_label,item]].dropna(how='any')
    figure(1)
    pl.title('Nightime ANN v Tower by year-month for '+item+' at '+Site_ID)

    try:
	xdata1a=Plottemp[item].groupby([lambda x: x.year,lambda x: x.month]).mean()
	plotxdata1a=True
    except:
	plotxdata1a=False
    try:
	xdata1b=Plottemp[NN_label].groupby([lambda x: x.year,lambda x: x.month]).mean()
	plotxdata1b=True
    except:
	plotxdata1b=False 
    if plotxdata1a==True:
	pl.plot(xdata1a,'r',label=item) 
    if plotxdata1b==True:
	pl.plot(xdata1b,'b',label=NN_label)
    pl.ylabel('Flux')    
    pl.xlabel('Year - Month')       
    pl.legend()
    pl.savefig(mypathforResults+'/ANN and Tower plots by year and month for variable '+item+' at '+Site_ID)
    #pl.show()
    pl.close()
    time.sleep(1)

def drawPr(tp,fp,tot,show=True):
    """
        draw the precision recall curve
    """
    det=numpy.array(sorted(tp+fp))
    atp=numpy.array(tp)
    afp=numpy.array(fp)
    #pylab.figure()
    #pylab.clf()
    rc=numpy.zeros(len(det))
    pr=numpy.zeros(len(det))
    #prc=0
    #ppr=1
    for i,p in enumerate(det):
        pr[i]=float(numpy.sum(atp>=p))/numpy.sum(det>=p)
        rc[i]=float(numpy.sum(atp>=p))/tot
        #print pr,rc,p
    ap=0
    for c in numpy.linspace(0,1,num=11):
        if len(pr[rc>=c])>0:
            p=numpy.max(pr[rc>=c])
        else:
            p=0
        ap=ap+p/11
    if show:
        pylab.plot(rc,pr,'-g')
        pylab.title("AP=%.3f"%(ap))
        pylab.xlabel("Recall")
        pylab.ylabel("Precision")
        pylab.grid()
        pylab.show()
        pylab.draw()
    return rc,pr,ap

def fixup_gauge(current_data):
    import pylab

    size = 34
    pylab.title("Pressure at Gauge 0", fontsize=size)
    pylab.xticks([1.0, 1.5, 2.0, 2.5, 3.0], fontsize=size)
    pylab.yticks([-0.3, -0.1, 0.1, 0.3, 0.5], fontsize=size)

def benchmark(clf_class, params, name):
    print("parameters:", params)
    t0 = time()
    clf = clf_class(**params).fit(X_train, y_train)
    print("done in %fs" % (time() - t0))

    if hasattr(clf, 'coef_'):
        print("Percentage of non zeros coef: %f"
              % (np.mean(clf.coef_ != 0) * 100))
    print("Predicting the outcomes of the testing set")
    t0 = time()
    pred = clf.predict(X_test)
    print("done in %fs" % (time() - t0))

    print("Classification report on test set for classifier:")
    print(clf)
    print()
    print(classification_report(y_test, pred,
                                target_names=news_test.target_names))

    cm = confusion_matrix(y_test, pred)
    print("Confusion matrix:")
    print(cm)

    # Show confusion matrix
    pl.matshow(cm)
    pl.title('Confusion matrix of the %s classifier' % name)
    pl.colorbar()

def makeContourPlot(scores, average, HEIGHT, WIDTH, outputId, maskId, plt_title, outputdir, barcodeId=-1, vmaxVal=100):
    pylab.bone()
    #majorFormatter = FormatStrFormatter('%.f %%')
    #ax = pylab.gca()
    #ax.xaxis.set_major_formatter(majorFormatter)
    
    pylab.figure()
    ax = pylab.gca()
    ax.set_xlabel(str(WIDTH) + ' wells')
    ax.set_ylabel(str(HEIGHT) + ' wells')
    ax.autoscale_view()
    pylab.jet()
    
    pylab.imshow(scores,vmin=0, vmax=vmaxVal, origin='lower')
    pylab.vmin = 0.0
    pylab.vmax = 100.0
    ticksVal = getTicksForMaxVal(vmaxVal)
    pylab.colorbar(format='%.0f %%',ticks=ticksVal)
    print "'%s'" % average
    if(barcodeId!=-1):
        if(barcodeId==0): maskId = "No Barcode Match,"
        else:             maskId = "Barcode Id %d," % barcodeId
    if plt_title != '': maskId = '%s\n%s' % (plt_title,maskId)
    print "Checkpoint A"
    pylab.title('%s Loading Density (Avg ~ %0.f%%)' % (maskId, average))
    pylab.axis('scaled')
    print "Checkpoint B"
    pngFn = outputdir+'/'+outputId+'_density_contour.png'
    print "Try save to", pngFn;
    pylab.savefig(pngFn, bbox_inches='tight')
    print "Plot saved to", pngFn;