def T2_cpmg_process(folder_to_process,plot='y'):
    """Given a folder of images will process cpmg data and return
    fitted T2 values and associated uncertainties"""
    data=img_roi_signal([folder_to_process],['EchoTime'])
    rois=data[0][0]
    TEs=data[2][0]
    mean_signal_mat=data[3]
    serr_signal_mat=data[4]
        
    if plot=='y':
        plt.figure()
    spin_echo_fits=[]
    for jj in np.arange(len(rois)-2):
        mean_sig=mean_signal_mat[0,jj,:]
        #serr_sig=serr_signal_mat[0,jj,:]
        mean_noise=np.mean(mean_signal_mat[0,-2,:])
        try:
            spin_echo_fit = SE_fit_new( np.array(TEs[0:]), mean_sig[0:], mean_noise, 'n' )
            if plot=='y':
                TE_full=np.arange(0,400,1)
                plt.subplot(4,4,jj+1)
                plt.plot(np.array(TEs[0:]), mean_sig[0:],'o')
                plt.plot(TE_full,spin_echo_fit(TE_full))
            
            spin_echo_fits.append(spin_echo_fit)
        except RuntimeError: 
            print 'RuntimeError'
            spin_echo=fitting.model('M0*exp(-x/T2)+a',{'M0':0,'T2':0,'a':0}) 
            spin_echo_fits.append(spin_echo)
    return spin_echo_fits   

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 embed_two_dimensions(data, vectorizer, size=10, n_components=5, colormap='YlOrRd'):
    if hasattr(data, '__iter__'):
        iterable = data
    else:
        raise Exception('ERROR: Input must be iterable')
    import itertools
    iterable_1, iterable_2 = itertools.tee(iterable)
    # get labels
    labels = []
    for graph in iterable_2:
        label = graph.graph.get('id', None)
        if label:
            labels.append(label)

    # transform iterable into sparse vectors
    data_matrix = vectorizer.transform(iterable_1)
    # embed high dimensional sparse vectors in 2D
    from sklearn import metrics
    distance_matrix = metrics.pairwise.pairwise_distances(data_matrix)

    from sklearn.manifold import MDS
    feature_map = MDS(n_components=n_components, dissimilarity='precomputed')
    explicit_data_matrix = feature_map.fit_transform(distance_matrix)

    from sklearn.decomposition import TruncatedSVD
    pca = TruncatedSVD(n_components=2)
    low_dimension_data_matrix = pca.fit_transform(explicit_data_matrix)

    plt.figure(figsize=(size, size))
    embed_dat_matrix_two_dimensions(low_dimension_data_matrix, labels=labels, density_colormap=colormap)
    plt.show()

def posterior(kpl, pk, err, pkfold=None, errfold=None):
    k0 = n.abs(kpl).argmin()
    kpl = kpl[k0:]
    if pkfold is None:
        print 'Folding for posterior'
        pkfold = pk[k0:].copy()
        errfold = err[k0:].copy()
        pkpos,errpos = pk[k0+1:].copy(), err[k0+1:].copy()
        pkneg,errneg = pk[k0-1:0:-1].copy(), err[k0-1:0:-1].copy()
        pkfold[1:] = (pkpos/errpos**2 + pkneg/errneg**2) / (1./errpos**2 + 1./errneg**2)
        errfold[1:] = n.sqrt(1./(1./errpos**2 + 1./errneg**2))

    #ind = n.logical_and(kpl>.2, kpl<.5)
    ind = n.logical_and(kpl>.15, kpl<.5)
    #ind = n.logical_and(kpl>.12, kpl<.5)
    #print kpl,pk.real,err
    kpl = kpl[ind]
    pk= kpl**3 * pkfold[ind]/(2*n.pi**2)
    err = kpl**3 * errfold[ind]/(2*n.pi**2)
    s = n.logspace(5.,6.5,100)
    data = []
    for ss in s:
        data.append(n.exp(-.5*n.sum((pk.real - ss)**2 / err**2)))
    #    print data[-1]
    data = n.array(data)
    #print data
    #print s
    #data/=n.sum(data)
    data /= n.max(data)
    p.figure(5)
    p.plot(s, data)
    p.plot(s, n.exp(-.5)*n.ones_like(s))
    p.plot(s, n.exp(-.5*2**2)*n.ones_like(s))
    p.show()

 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 trace(data, name, format='png', datarange=(None, None), suffix='', path='./', rows=1, columns=1, 
    num=1, last=True, fontmap = None, verbose=1):
    """
    Generates trace plot from an array of data.

    :Arguments:
        data: array or list
            Usually a trace from an MCMC sample.

        name: string
            The name of the trace.
            
        datarange: tuple or list
            Preferred y-range of trace (defaults to (None,None)).

        format (optional): string
            Graphic output format (defaults to png).

        suffix (optional): string
            Filename suffix.

        path (optional): string
            Specifies location for saving plots (defaults to local directory).
            
        fontmap (optional): dict
            Font map for plot.

    """

    if fontmap is None: fontmap = {1:10, 2:8, 3:6, 4:5, 5:4}

    # Stand-alone plot or subplot?
    standalone = rows==1 and columns==1 and num==1

    if standalone:
        if verbose>0:
            print_('Plotting', name)
        figure()

    subplot(rows, columns, num)
    pyplot(data.tolist())
    ylim(datarange)

    # Plot options
    title('\n\n   %s trace'%name, x=0., y=1., ha='left', va='top', fontsize='small')

    # Smaller tick labels
    tlabels = gca().get_xticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    tlabels = gca().get_yticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    if standalone:
        if not os.path.exists(path):
            os.mkdir(path)
        if not path.endswith('/'):
            path += '/'
        # Save to file
        savefig("%s%s%s.%s" % (path, name, suffix, format))

    def plot_matches(self, name, show_below = True, match_maximum = None):
        """ 対応点を線で結んで画像を表示する
          入力： im1,im2（配列形式の画像）、locs1,locs2（特徴点座標）
             machescores（match()の出力）、
             show_below（対応の下に画像を表示するならTrue）"""
        im1 = self._image_1.get_array_image()
        im2 = self._image_2.get_array_image()
        self.appendimages()
        im3 = self._append_image
        if self._match_score is None:
            self.match()
        locs1 = self._image_1.get_shift_location()
        locs2 = self._image_2.get_shift_location()
        if show_below:
            im3 = numpy.vstack((im3,im3))
        pylab.figure(dpi=160)
        pylab.gray()
        pylab.imshow(im3, aspect = 'auto')

        cols1 = im1.shape[1]
        match_num = 0
        for i,m in enumerate(self._match_score):
            if m > 0 : 
                pylab.plot([locs1[i][0],locs2[m][0]+cols1], [locs1[i][1],locs2[m][1]], 'c')
                match_num = match_num + 1
            if match_maximum is not None and match_num >= match_maximum:
                break
        pylab.axis('off')
        pylab.savefig(name, dpi=160)

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 cmap_plot(cmdLine):

    pylab.figure(figsize=[5,10])
    a=outer(ones(10),arange(0,1,0.01))
    subplots_adjust(top=0.99,bottom=0.00,left=0.01,right=0.8)
    maps=[m for m in cm.datad if not m.endswith("_r")]
    maps.sort()
    l=len(maps)+1
    for i, m in enumerate(maps):
        print m
        subplot(l,1,i+1)
        pylab.setp(pylab.gca(),xticklabels=[],xticks=[],yticklabels=[],yticks=[])
        imshow(a,aspect='auto',cmap=get_cmap(m),origin="lower")
        pylab.text(100.85,0.5,m,fontsize=10)

# render plot

    if cmdLine: 
        pylab.show(block=True)
    else: 
        pylab.ion()
        pylab.plot([])
        pylab.ioff()
	
    status = 1
    return status

def plotAllWarmJumps():
    jumpAddrs = np.array(getAllWarmJumpsAddr()).reshape((8, 18))
    figure()
    pcolor(jumpAddrs)
    for (x, y), v in np.ndenumerate(jumpAddrs):
        text(y + 0.125, x + 0.5, "0x%03x" % v)
    show()

    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 plot_heatingrate(data_dict, filename, do_show=True):
    pl.figure(201)
    color_list = ['b','r','g','k','y','r','g','b','k','y','r',]
    fmtlist = ['s','d','o','s','d','o','s','d','o','s','d','o']
    result_dict = {}
    for key in data_dict.keys():
        x = data_dict[key][0]
        y = data_dict[key][1][:,0]
        y_err = data_dict[key][1][:,1]

        p0 = np.polyfit(x,y,1)
        fit = LinFit(np.array([x,y,y_err]).transpose(), show_graph=False)
        p1 = [0,0]
        p1[0] = fit.param_dict[0]['Slope'][0]
        p1[1] = fit.param_dict[0]['Offset'][0]
        print fit
        x0 = np.linspace(0,max(x))
        cstr = color_list.pop(0)
        fstr = fmtlist.pop(0)
        lstr = key + " heating: {0:.2f} ph/ms".format((p1[0]*1e3)) 
        pl.errorbar(x/1e3,y,y_err,fmt=fstr + cstr,label=lstr)
        pl.plot(x0/1e3,np.polyval(p0,x0),cstr)
        pl.plot(x0/1e3,np.polyval(p1,x0),cstr)
        result_dict[key] = 1e3*np.array(fit.param_dict[0]['Slope'])
    pl.xlabel('Heating time (ms)')
    pl.ylabel('nbar')
    if do_show:
        pl.legend()
        pl.show()
    if filename != None:
        pl.savefig(filename)
    return result_dict

    def swi_histogram(self,dir,name,measure,dpi=80,width=8,height=6,b_left='0.1',b_bot='0.1',b_top='0.1',b_right='0.1',bins='20'):
        s=ccm.stats.Stats('%s/%s'%(dir,name))
        data=s.get_raw(measure) 
        
        bins=int(bins)
        
           

        pylab.figure(figsize=(float(width),float(height)))
        try: b_left=float(b_left)
        except: b_left=0.1
        try: b_right=float(b_right)
        except: b_right=0.1
        try: b_top=float(b_top)
        except: b_top=0.1
        try: b_bot=float(b_bot)
        except: b_bot=0.1
        pylab.axes((b_left,b_bot,1.0-b_left-b_right,1.0-b_top-b_bot))
        
        
        pylab.hist(data,bins=bins)

            
        
        img=StringIO.StringIO()
        if type(dpi) is list: dpi=dpi[-1]
        pylab.savefig(img,dpi=int(dpi),format='png')
        return 'image/png',img.getvalue()

def T1_ir_bootstrap(folders_to_process,N=1000,plot='n'):    
    """Given a folder of images will process IR data and return
    fitted T1 values and associated uncertainties.  Uncertainties 
    are obtained by bootstrapping pixels within each ROI"""  
    images=[read_dicoms(folder,['InversionTime']) for folder in folders_to_process]
    all_rois=[load_ROIs(folder+'/rois') for folder in folders_to_process]   
    TI_images=[img[0][0] for img in images]
    TIs=np.array([img[1][0]['InversionTime'] for img in images])
    T1s=[]
    T1_errs=[]
    for kk in np.arange(len(all_rois[0])-2):
        rois=[roi_list[kk] for roi_list in all_rois]
        sig=np.zeros(len(rois))
        T1bs=[]
        for mm, roi in enumerate(rois):
            sig[mm]=TI_images[mm][roi.get_indices()].mean()
        fit = IR_fit(TIs, sig) 
        T1s.append(fit['T1'].value)
        if plot=='y':
            plt.figure()
            plt.plot(TIs,sig,'o')
            TI_full=np.arange(0,8000,10)
            plt.plot(TI_full,fit(TI_full))               
        if N>0:  
            for nn in np.arange(N):
                for mm,roi in enumerate(rois):
                    npix=len(roi.get_indices()[0])
                    pixels=roi.get_indices()
                    ind=np.random.randint(npix,size=npix)
                    sig[mm]=TI_images[mm][pixels[0][ind],pixels[1][ind]].mean()                
                fit = IR_fit(TIs, sig) 
                T1bs.append(fit['T1'].value)
            #T1s.append(np.mean(T1bs))
            T1_errs.append(np.std(T1bs))
    return T1s, T1_errs

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 plot_sphere_x( s, fname ):
  """ put plot of ionization fractions from sphere `s` into fname """

  plt.figure()
  s.Edges.units = 'kpc'
  s.r_c.units = 'kpc'
  xx = s.r_c
  L = s.Edges[-1]

  plt.plot( xx, np.log10( s.xHe1 ),
            color='green', ls='-', label = r'$x_{\rm HeI}$' )
  plt.plot( xx, np.log10( s.xHe2 ),
            color='green', ls='--', label = r'$x_{\rm HeII}$' )
  plt.plot( xx, np.log10( s.xHe3 ),
            color='green', ls=':', label = r'$x_{\rm HeIII}$' )

  plt.plot( xx, np.log10( s.xH1 ),
            color='red', ls='-', label = r'$x_{\rm HI}$' )
  plt.plot( xx, np.log10( s.xH2 ),
            color='red', ls='--', label = r'$x_{\rm HII}$' )

  plt.xlim( -L/20, L+L/20 )
  plt.xlabel( 'r_c [kpc]' )

  plt.ylim( -4.5, 0.2 )
  plt.ylabel( 'log 10 ( x )' )

  plt.grid()
  plt.legend(loc='best', ncol=2)
  plt.tight_layout()
  plt.savefig( 'doc/img/x_' + fname )

    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 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;

def geweke_plot(data, name, format='png', suffix='-diagnostic', path='./', fontmap = None, 
    verbose=1):
    # Generate Geweke (1992) diagnostic plots

    if fontmap is None: fontmap = {1:10, 2:8, 3:6, 4:5, 5:4}

    # Generate new scatter plot
    figure()
    x, y = transpose(data)
    scatter(x.tolist(), y.tolist())

    # Plot options
    xlabel('First iteration', fontsize='x-small')
    ylabel('Z-score for %s' % name, fontsize='x-small')

    # Plot lines at +/- 2 sd from zero
    pyplot((nmin(x), nmax(x)), (2, 2), '--')
    pyplot((nmin(x), nmax(x)), (-2, -2), '--')

    # Set plot bound
    ylim(min(-2.5, nmin(y)), max(2.5, nmax(y)))
    xlim(0, nmax(x))

    # Save to file
    if not os.path.exists(path):
        os.mkdir(path)
    if not path.endswith('/'):
        path += '/'
    savefig("%s%s%s.%s" % (path, name, suffix, format))

    def displayResults(self,res, cm=pylab.cm.gray, title='Specify a title'):
        if self.display:
		self.count=self.count+1
        	pylab.figure(self.count)
        	pylab.imshow(res, cm, interpolation='nearest')
        	pylab.colorbar()
        	pylab.title(title)