def scatter_stats(db, s1, s2, f1=None, f2=None, **kwargs):
    if f1 == None:
        f1 = lambda x: x  # constant function

    if f2 == None:
        f2 = f1

    x = []
    xerr = []

    y = []
    yerr = []

    for k in db:
        x_k = [f1(x_ki) for x_ki in db[k].__getattribute__(s1).gettrace()]
        y_k = [f2(y_ki) for y_ki in db[k].__getattribute__(s2).gettrace()]

        x.append(pl.mean(x_k))
        xerr.append(pl.std(x_k))

        y.append(pl.mean(y_k))
        yerr.append(pl.std(y_k))

        pl.text(x[-1], y[-1], " %s" % k, fontsize=8, alpha=0.4, zorder=-1)

    default_args = {"fmt": "o", "ms": 10}
    default_args.update(kwargs)
    pl.errorbar(x, y, xerr=xerr, yerr=yerr, **default_args)
    pl.xlabel(s1)
    pl.ylabel(s2)

def data_to_ch(data):
    ch = {}
    for ch_ind in range(1, 97):
        ch[ch_ind] = {}
        ch[ch_ind]["bl"] = data[ch_ind]["blanks"]
        ch[ch_ind]["bl_mu"] = pl.mean(ch[ch_ind]["bl"])
        ch[ch_ind]["bl_sem"] = pl.std(ch[ch_ind]["bl"]) / pl.sqrt(len(ch[ch_ind]["bl"]))
        for ind in sorted(data[ch_ind].keys()):
            if ind != "blanks":
                k = ind[0]
                if k not in ch[ch_ind]:
                    ch[ch_ind][k] = {}
                    ch[ch_ind][k]["fr"] = []
                    ch[ch_ind][k]["fr_mu"] = []
                    ch[ch_ind][k]["fr_sem"] = []
                    ch[ch_ind][k]["pos_y"] = []
                    ch[ch_ind][k]["dprime"] = []
                ch[ch_ind][k]["fr"].append(data[ch_ind][ind]["on"])
                ch[ch_ind][k]["fr_mu"].append(pl.mean(data[ch_ind][ind]["on"]))
                ch[ch_ind][k]["fr_sem"].append(pl.std(data[ch_ind][ind]["on"]) / pl.sqrt(len(data[1][ind]["on"])))
                ch[ch_ind][k]["pos_y"].append(ind[2])
                # print ch[ch_ind][k]['pos_y']
                # print pl.std(data[ch_ind][ind]['on'])
                ch[ch_ind][k]["dprime"].append(
                    (pl.mean(data[ch_ind][ind]["on"]) - ch[ch_ind]["bl_mu"])
                    / ((pl.std(ch[ch_ind]["bl"]) + pl.std(data[ch_ind][ind]["on"])) / 2)
                )
                # print ch[ch_ind]['OSImage_5']['pos_y']
    return ch

def compare_models(db, stoch="itn coverage", stat_func=None, plot_type="", **kwargs):
    if stat_func == None:
        stat_func = lambda x: x

    X = {}
    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c] = []

    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c].append([stat_func(x_ki) for x_ki in db[k].__getattribute__(stoch).gettrace()])

    x = pl.array([pl.mean(xc[0]) for xc in X.values()])
    xerr = pl.array([pl.std(xc[0]) for xc in X.values()])
    y = pl.array([pl.mean(xc[1]) for xc in X.values()])
    yerr = pl.array([pl.std(xc[1]) for xc in X.values()])

    if plot_type == "scatter":
        default_args = {"fmt": "o", "ms": 10}
        default_args.update(kwargs)
        for c in X.keys():
            pl.text(pl.mean(X[c][0]), pl.mean(X[c][1]), " %s" % c, fontsize=8, alpha=0.4, zorder=-1)
        pl.errorbar(x, y, xerr=xerr, yerr=yerr, **default_args)
        pl.xlabel("First Model")
        pl.ylabel("Second Model")
        pl.plot([0, 1], [0, 1], alpha=0.5, linestyle="--", color="k", linewidth=2)

    elif plot_type == "rel_diff":
        d1 = sorted(100 * (x - y) / x)
        d2 = sorted(100 * (xerr - yerr) / xerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean (%)")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err (%)")
    elif plot_type == "abs_diff":
        d1 = sorted(x - y)
        d2 = sorted(xerr - yerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err")
    else:
        assert 0, "plot_type must be abs_diff, rel_diff, or scatter"

    return pl.array([x, y, xerr, yerr])

def plot2():
    import pylab as pl
    hs, ds = [], []
    for event, time in load():
        if event == main_start:
            start_time = time
        elif event == main_end:
            d0, h0 = days_hours(start_time)
            d1, h1 = days_hours(time)
            hs.append((h0, h1))
            ds.append((d0, d1))
            pl.plot([d0, d1], [h0, h1], 'b')
    ihs, fhs = zip(*hs)
    ids, fds = zip(*ds)
    pl.plot(ids, ihs, 'g')
    pl.plot([ids[0], ids[-1]], [pl.mean(ihs)] * 2, 'g--')
    pl.plot(fds, fhs, 'r')
    pl.plot([fds[0], fds[-1]], [pl.mean(fhs)] * 2, 'r--')
    f, i = pl.mean(fhs), pl.mean(ihs)
    pl.plot([fds[0], fds[-1]], [(f + i) / 2] * 2, 'b--')
    print i, f, f - i, (f + i) / 2
    std_i, std_f = pl.std(ihs), pl.std(fhs)
    print std_i, std_f
    pl.xlim(ids[0], fds[-1])
    pl.ylim(4, 28)
    pl.grid(True)
    pl.xlabel('Time [day]')
    pl.ylabel('Day interval [hours]')
    pl.show()

 def stderr(X,Y=None):
     if len(X) <= 1: return 0.0
     stderr_x = pow(pylab.std(X),2)/len(X)
     if Y:
         if len(Y) <= 1: return 0.0
         stderr_y = pow(pylab.std(Y),2)/len(Y)
     else: stderr_y = 0
     return math.sqrt(stderr_x + stderr_y)

def Vplot(Ws):
	"""Calculate the potential function and plot it"""

	N_bstrp = input("Please enter the number of bootstraps: ")
	N_bin = input("Please enter the bin size: ")
	style = raw_input("Please enter a linestyle: ")

	Ws = bin(Ws,N_bin)
	aVs = pl.zeros((N_bstrp,) + pl.shape(Ws)[1:])
	bs = pl.zeros((N_bstrp,3))
        
	for i in xrange(N_bstrp):
		W = pl.mean(bootstrap(Ws),axis=0)
		aVs[i] = calcaV(W,method="fit")
		bs[i] = potfit(aVs[i,:,0])

			
	r = pl.arange(1,7)
	aV = pl.mean(aVs,axis=0)
	aVerr = pl.std(aVs,axis=0)
	b = pl.mean(bs,axis=0)

	a_s = 0.5 / pl.sqrt((1.65 + bs[:,1]) / bs[:,0])
	sigmas = bs[:,0] / a_s**2
	Bs = bs[:,1]
	As = bs[:,2] / a_s

	a = pl.mean(a_s)
	aerr = pl.std(a_s)
	sigma = pl.mean(sigmas)
	sigmaerr = pl.std(sigmas)
	B = pl.mean(Bs)
	Berr = pl.std(Bs)
	A = pl.mean(As)
	Aerr = pl.std(As)

	print("Fit parameters:")
	print("sigma = %f +/- %f fm^-2 = %f +/- %f MeV^2"
		% (sigma, sigmaerr, sigma * 197**2, sigmaerr * 197**2))
	print("B = %f +/- %f" % (B, Berr))
	print("A = %f +/- %f fm^-1 = %f +/- %f MeV"
		% (A, Aerr, A*197, Aerr*197))
	print("Lattice spacing, a = %f +/- %f fm = %f +/- %f MeV^-1"
		% (a, aerr, a/197, aerr/197))
	
	r_fit = pl.arange(0.25,r[-1]+1,0.1)
	aV_fit = V(b,r_fit)
	
	handles = []
	handles.append(pl.errorbar(r,aV[:,0],yerr=aVerr[:,0],fmt='o'+style[0]))
	handles.append(pl.plot(r_fit,aV_fit,style))
	pl.ylim([0,pl.nanmax(aV)+0.25])
	pl.xlim([0,pl.nanmax(r_fit)+0.25])
	pl.xlabel("$r / a$")
	pl.ylabel("$aV(r)$")

	return aV,handles

 def update(self,t,val):
     oldavg = self.avg
     avrg = Ema.update(self,t,val)
     
     self.__samples.append((t,self.lastvalue))
     self.__samples_nonone.append(self.lastvalue)
     
     if oldavg == None:
         self.first_t = t 
         return (None,None,None)
     newavg = avrg
     
     # check limits of timeframe 
     while t - self.__samples[0][0] > self.timeframe and len(self.__samples) > 2:
         _, pv = self.__samples.pop(0)
         del self.__samples_nonone[0]
         #self.variance += (val-pv)*(val-newavg+pv-oldavg)/(self.timeframe) # this seems to use constant number of samples
         #std = math.sqrt(self.variance)
         
     std = pylab.std(self.__samples_nonone)*self.k # this takes long
     
     if avrg != None:
         return avrg,avrg+std,avrg-std
     else:
         return (None,None,None)

def flow_rate_hist(sheets):
    ant_rates = []
    weights = []
    for sheet in sheets:
        ants, seconds, weight = flow_rate(sheet)
        ant_rate = seconds / ants
        #ant_rate = ants / seconds
        ant_rates.append(ant_rate)
        weights.append(float(weight))
        #weights.append(seconds)

    weights = pylab.array(weights)
    weights /= sum(weights)

    #print "ants per second"
    print "seconds per ant"
    mu = pylab.mean(ant_rates)
    print "mean", pylab.mean(ant_rates)
    wmean = pylab.average(ant_rates, weights=weights)
    print "weighted mean", wmean
    print "median", pylab.median(ant_rates)
    print "std", pylab.std(ant_rates, ddof=1)
    ant_rates = pylab.array(ant_rates)
    werror = (ant_rates - mu) * weights
    print "weighted std", ((sum(werror ** 2))) ** 0.5
    print "weighted std 2", (pylab.average((ant_rates - mu)**2, weights=weights)) ** 0.5
    pylab.figure()
    pylab.hist(ant_rates)
    pylab.savefig('ant_flow_rates.pdf', format='pdf')
    pylab.close()

def _CalcMutualNearestNeighbors(hull_points, all_points):
    all_points_list = list(all_points)
    ds = distance.pdist(list(all_points))
    std_d = p.std(ds)
    
    square_ds = distance.squareform(ds)
    nearest_neighbors = {}
    
    for i, point in enumerate(all_points_list):
        if point not in hull_points:
            continue
        
        my_ds = [(d, j) for j, d in enumerate(square_ds[i])
                 if j != i]
        my_ds.sort()
        nearest_neighbors[point] = set([j for d,j in my_ds[:3]])
    
    no_mutual = set()
    for i, point in enumerate(all_points_list):
        if point not in hull_points:
            continue
        
        no_nbrs = True
        for neighbor_index in nearest_neighbors.get(point, []):
            neighbor = all_points_list[neighbor_index]
            neighbor_set = nearest_neighbors.get(neighbor, [])
            if i in neighbor_set:
                no_nbrs = False
        
        if no_nbrs:
            no_mutual.add(point)
                
    return no_mutual

 def generate_normalized_test_data(self,
                                   channels, 
                                   time_points, 
                                   function, 
                                   sampling_frequency, 
                                   initial_phase=0.0):
     """
     A method which generates a normalized (mu = 0, sigma =1) signal for testing, with
     the specified number of "channels" which are all generated using the given function
     """
    
     #Generate an empty ndarray
     data = numpy.zeros((time_points, channels))
     
     #Compute the values for all channels
     for channel_index in range(channels):
         for time_index in range(time_points):
             data[time_index, channel_index] = function(2.0 * numpy.pi * (channel_index + 1) * (time_index / sampling_frequency + initial_phase))
         current_channel = data[:, channel_index]
         current_channel = (current_channel - pylab.mean(current_channel))/pylab.std(current_channel)
         data[:, channel_index] = current_channel
         
     #Generate a time series build out of the data
     test_data = TimeSeries(input_array = data, 
                            channel_names = [("test_channel_%s" % i) for i in range(channels)],
                            sampling_frequency = sampling_frequency,
                            start_time = initial_phase,
                            end_time = float(time_points) / sampling_frequency + initial_phase)
     
     return test_data

def pc_pm_std(data, ndim):
    """
    This is a helper function.
    It returns the value of +1 * std(x), where x is the ndim-th principal
    component of the data

    Parameters:
    -----------
    data: `array` (*n*-by-*d*)
        the data on which the principal component analysis is performed.
    ndim: `integer`
        the number of the principal axis on which the analysis is performed.
        **NOTE** this is zero-based, i.e. to compute the first principal
        component, ndim=0

    Returns:
    --------
    std_pc: `array` (1-by-*d*)
        the vector that points in the direction of the *ndim*th principal
        axis, and has the length of the standard deviation of the scores
        along this axis.

    """

    u,s,v = svd(data.T, full_matrices = False)
    direction = u[:, ndim : ndim + 1]
    scale = std(dot(direction.T, data.T))
    return scale * direction.T

    def zeroPaddData(self,desiredLength,paddmode='zero',where='end'):    
        #zero padds the time domain data, it is possible to padd at the beginning,
        #or at the end, and further gaussian or real zero padding is possible        
        #might not work for gaussian mode!

        desiredLength=int(desiredLength)
        #escape the function        
        if desiredLength<0:
            return 0

        #calculate the paddvectors        
        if paddmode=='gaussian':
            paddvec=py.normal(0,py.std(self.getPreceedingNoise())*0.05,desiredLength)
        else:
            paddvec=py.ones((desiredLength,self.tdData.shape[1]-1))
            paddvec*=py.mean(self.tdData[-20:,1:])
            
        timevec=self.getTimes()
        if where=='end':
            #timeaxis:
            newtimes=py.linspace(timevec[-1],timevec[-1]+desiredLength*self.dt,desiredLength)
            paddvec=py.column_stack((newtimes,paddvec))
            longvec=py.row_stack((self.tdData,paddvec))
        else:
            newtimes=py.linspace(timevec[0]-(desiredLength+1)*self.dt,timevec[0],desiredLength)
            paddvec=py.column_stack((newtimes,paddvec))
            longvec=py.row_stack((paddvec,self.tdData))
            
        self.setTDData(longvec)

 def post_lecture(self):
     STD = std(self.Y, 1)
     MM = mean(self.Y, 1)
     TT, self.NN = self.Y.shape
     if self.centred:
         for t in xrange(0, TT):
             self.Y[t, :] = (self.Y[t, :] - MM[t]) / STD[t]

 def getelnNoise(self,tdData):
     #returns the uncertainty due to electronic noise
     
     #signal preceeding the pulse (X and Y channel)  
     precNoise=self.getPreceedingNoise(tdData)
     #is this normalization really correct?!
     elnNoise = py.std(precNoise, ddof = 1,axis=0)/py.sqrt(precNoise.shape[0])
     return elnNoise

def readDatDirectory(key, directory):
    global stats
    #Don't read data in if it's already read
    if not key in DATA["mean"]:
        data = defaultdict(array)

        #Process the dat files
        for datfile in glob.glob(directory + "/*.dat"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Process the div files'
        for datfile in glob.glob(directory + "/*.div"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Iterate through the stats and calculate mean/standard deviation
        for aKey in stats:
            if aKey in data:
                DATA["mean"][key][aKey] = mean(data[aKey], axis=0)
                DATA["median"][key][aKey] = median(data[aKey], axis=0)
                DATA["std"][key][aKey] = std(data[aKey], axis=0)
                DATA["ste"][key][aKey] = std(data[aKey], axis=0)/ sqrt(len(data[aKey]))
                DATA["min"][key][aKey] = mean(data[aKey], axis=0)-amin(data[aKey], axis=0)
                DATA["max"][key][aKey] = amax(data[aKey], axis=0)-mean(data[aKey], axis=0)
                DATA["actual"][key][aKey] = data[aKey]

def xyamb(xytab,qu,xyout=''):

    mytb=taskinit.tbtool()

    if not isinstance(qu,tuple):
        raise Exception,'qu must be a tuple: (Q,U)'

    if xyout=='':
        xyout=xytab
    if xyout!=xytab:
        os.system('cp -r '+xytab+' '+xyout)

    QUexp=complex(qu[0],qu[1])
    print 'Expected QU = ',qu   # , '  (',pl.angle(QUexp)*180/pi,')'

    mytb.open(xyout,nomodify=False)

    QU=mytb.getkeyword('QU')['QU']
    P=pl.sqrt(QU[0,:]**2+QU[1,:]**2)

    nspw=P.shape[0]
    for ispw in range(nspw):
        st=mytb.query('SPECTRAL_WINDOW_ID=='+str(ispw))
        if (st.nrows()>0):
            q=QU[0,ispw]
            u=QU[1,ispw]
            qufound=complex(q,u)
            c=st.getcol('CPARAM')
            fl=st.getcol('FLAG')
            xyph0=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
            print 'Spw = '+str(ispw)+': Found QU = '+str(QU[:,ispw])  # +'   ('+str(pl.angle(qufound)*180/pi)+')'
            #if ( (abs(q)>0.0 and abs(qu[0])>0.0 and (q/qu[0])<0.0) or
            #     (abs(u)>0.0 and abs(qu[1])>0.0 and (u/qu[1])<0.0) ):
            if ( pl.absolute(pl.angle(qufound/QUexp)*180/pi)>90.0 ):
                c[0,:,:]*=-1.0
                xyph1=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
                st.putcol('CPARAM',c)
                QU[:,ispw]*=-1
                print '   ...CONVERTING X-Y phase from '+str(xyph0)+' to '+str(xyph1)+' deg'
            else:
                print '      ...KEEPING X-Y phase '+str(xyph0)+' deg'
            st.close()
    QUr={}
    QUr['QU']=QU
    mytb.putkeyword('QU',QUr)
    mytb.close()
    QUm=pl.mean(QU[:,P>0],1)
    QUe=pl.std(QU[:,P>0],1)
    Pm=pl.sqrt(QUm[0]**2+QUm[1]**2)
    Xm=0.5*atan2(QUm[1],QUm[0])*180/pi

    print 'Ambiguity resolved (spw mean): Q=',QUm[0],'U=',QUm[1],'(rms=',QUe[0],QUe[1],')','P=',Pm,'X=',Xm

    stokes=[1.0,QUm[0],QUm[1],0.0]
    print 'Returning the following Stokes vector: '+str(stokes)
    
    return stokes

  def __init__(self, fp, delimiter="\t", require_header=False):
    """Load matrix of floats into self.

    File Format:
    ===========
    First line: Column titles. 
    -------
      First column is row variable name, all other columns are sample IDs
      e.g.:
        miRNA_ID	sample_id_1	sample_id_2...
        
    Next lines: rows of float data. Each row represents a single variable.
    -----
      e.g. 
        hsa-let-7e	332.0	690.0...

    Args:
      fp: [*str] like open filepointer to matrix of floats.
      delimiter: str of column delimiter
      require_header: bool if to require first line file header from input
    """
    self.rows = {}

    for line in fp:
      # Skip header comments
      if line[0] == "#": continue
      first_row = line.strip('\n').split(delimiter)
      
      # Verify that the first line looks like column headers.
      if require_header and first_row != "miRNA_ID":
        # Headers are required but the first row seems malformed. Error.
        raise ValueError, "Line not valid matrix column table header."
      elif first_row[0] == "miRNA_ID":
        # This looks a list of sample_IDs as column titles. Set samples.
        self.samples = first_row[1:]
      else:
        # This looks like a data row. Add it.
        self.samples = None
        self._add_row(first_row)
      break # Exit loop after parsing first non-comment line

    # Parse all other lines as rows of floats named by the first column entry.
    for line in fp:
      # WARNING: DO NOT STRIP TRAILING TABS
      row = line.strip('\n').split(delimiter)
      self._add_row(row)

    # Set matrix dimensions from last row.
    self.n = len(row) -1 # first entry is variable name
    self.m = len(self.rows)

    # Compute and store all variable standard deviations.
    self.stds = {}
    for name, values in self.rows.items():
      # Remove None's from values; DO NOT remove zeros!
      std = pylab.std(filter(lambda x: x is not None, values)) 
      self.stds[name] = std

def plot_histogram(histogram, html_writer, title='', max_pathway_length=8, xmin=None, xlim=20, error_bars=True, min_to_show=20, legend_loc='upper left'):
    fig = pylab.figure()

    pylab.hold(True)

    reps = 1000
    
    y_offset = 0
    offset_step = 0.007
    colors = {1:'r', 2:'orange', 3:'green', 4:'cyan', 5:'blue', 'Rest':'violet', 'Not first':'k--', 'No known regulation':'grey', 'Activated':'green', 'Inhibited':'r', 'Mixed regulation':'blue'}
    for key, value in histogram.iteritems():
        if len(value) >= min_to_show:
            m = stats.cmedian(value)
            
            sample_std = None
            
            if error_bars:
                sample_vals = []
                i = 0
                while i < reps:
                    samples = []
                    while len(samples) < len(value):
                        samples.append(random.choice(value))
                    sample_vals.append(pylab.median(samples))
                    i += 1
                
                sample_std = pylab.std(sample_vals)
                        
            plotting.cdf(value, label='%s (med=%.1f, N=%d)' % \
                (key, m, len(value)),
                style=colors.get(key, 'grey'), std=sample_std, y_offset=y_offset)
            y_offset += offset_step
            

    xmin = -1 * xlim if xmin == None else xmin
    pylab.xlim(xmin, xlim)
    pylab.xlabel('Irreversability')
    #pylab.xlabel('deltaG')
    pylab.ylabel('Cumulative distribution')
    legendfont = matplotlib.font_manager.FontProperties(size=11)
    pylab.legend(loc=legend_loc, prop=legendfont)
    pylab.title(title)
    pylab.hold(False)
    
    if 'Not first' in histogram:
        print '%s, first vs. non-first ranksum test: ' % title + '(%f, %f)' % stats.ranksums(histogram[1], histogram['Not first'])
    
    if 'Inhibited' in histogram:
        print '%s, inhibited vs. non-regulated ranksum test: ' % title + '(%f, %f)' % stats.ranksums(histogram['Inhibited'], histogram['No known regulation'])
         
    
    #for k1, h1 in histogram.iteritems():
    #    for k2, h2 in histogram.iteritems():
    #        print k1, k2, stats.ranksums(h1, h2)
    
    return fig

    def old_bar_graph_WOULD_NEED_FIX(
           img_dst, points, independent, dependent, search):
        from pylab import arange, xticks, xlim, mean, std
        # Lines
        plots = []
        labels = []
        args = {}
        for point in points:
            args[point.get(independent)] = True
        keys = args.keys()
        keys.sort()
        ind = {}
        bars = {}
        for k in keys:
            ind[k] = keys.index(k)
            bars[k] = {}
        # Float is necessary
        ticks = arange(len(args), dtype='float')

        for point_type in selected:
            x_val = []
            y_val = []
            for point in selected[point_type]:
                x_val.append(point.get(independent))
                y_val.append(point.get(dependent))

            width = 1./(len(selected)+2)
            for point in selected[point_type]:
                pi = point.get(independent)
                pd = point.get(dependent)
                # bars[independent][screen] = ([points],position,color)
                if point_type not in bars[pi].keys():
                    bars[pi][point_type] = ([pd], \
                                            ind[pi]+width*len(bars[pi]), \
                                            color)
                else:
                    bars[pi][point_type][0].append(pd)
            labels.append(label_type(point_type))

        plotlist = {}
        for pi in bars:
            for pt in bars[pi]:
                if len(bars[pi][pt][0]) > 1:
                    error = std(bars[pi][pt][0])
                else:
                    error = 0
                p = bar(bars[pi][pt][1],mean(bars[pi][pt][0]),width,\
                      color=bars[pi][pt][2], yerr=error)
                plotlist[pt] = p
                ticks[ind[pi]] += width/2
        plots = plotlist.values()
        plots.sort()
        keys = args.keys()
        keys.sort()
        xticks(ticks, keys)
        xlim(-width,len(ticks))

 def __init__(self, v, bin_precision=5):
   self.size = len(v)
   self.mean = pylab.mean(v)
   self.std = pylab.std(v)
   self.min = min(v)
   self.max = min(v)
   self.bins = {}
   for x in v:
     key = ("%%.%df" % precision) % x 
     bins[key] = bins.get(key, 0) + 1