def fade( from_color, to_color, steps ):
    from_color = narray(from_color)
    to_color = narray(to_color)
    diff = (to_color-from_color)/float(steps)
    for n in range(steps):
        color = from_color + n*diff
        yield color.round().astype(int).tolist()

 def read(self,file):
     iclose=False
     try:
         # If file is not open, open it as a fileobject
         fo=open(file,'rb')
         iclose=True
     except:
         # If already a fileobject, we will append to it
         fo=file
     ti=array('i')
     ti.read(fo,2)
     self.nboot=ti[0]
     self.ndim=ti[1]
     self.values.resize(self.nboot,self.ndim)
     tf=array('d')
     tf.read(fo,self.ndim)
     self.Avg=narray(tf)
     tf=array('d')
     tf.read(fo,self.ndim)
     self.Std=narray(tf)
     for iboot in xrange(self.nboot):
         tf=array('d')
         tf.read(fo,self.ndim)
         self.values[iboot]=narray(tf)
     
     # If file was a new file, close it
     if iclose:
         fo.close()

def _blob(x,y,area,colour):
	"""
	Draws a square-shaped blob with the given area (< 1) at
	the given coordinates.
	"""
	hs = sqrt(area) / 2
	xcorners = narray([x - hs, x + hs, x + hs, x - hs])
	ycorners = narray([y - hs, y - hs, y + hs, y + hs])
	P.fill(xcorners, ycorners, colour, edgecolor=colour)

    def tolist(self, fill_value=None):
        """
        Return the data portion of the array as a list.

        Data items are converted to the nearest compatible Python type.
        Masked values are converted to fill_value. If fill_value is None,
        the corresponding entries in the output list will be ``None``.

        """
        if fill_value is not None:
            return self.filled(fill_value).tolist()
        result = narray(self.filled().tolist(), dtype=object)
        mask = narray(self._mask.tolist())
        result[mask] = None
        return result.tolist()

 def __array_finalize__(self,obj):
     # Make sure we have a _fieldmask by default ..
     _fieldmask = getattr(obj, '_fieldmask', None)
     if _fieldmask is None:
         mdescr = _make_mask_dtype(ndarray.__getattribute__(self, 'dtype'))
         _mask = getattr(obj, '_mask', nomask)
         if _mask is nomask:
             _fieldmask = np.empty(self.shape, dtype=mdescr).view(recarray)
             _fieldmask.flat = tuple([False]*len(mdescr))
         else:
             _fieldmask = narray([tuple([m]*len(mdescr)) for m in _mask],
                                 dtype=mdescr).view(recarray)
     # Update some of the attributes
     if obj is not None:
         _baseclass = getattr(obj,'_baseclass',type(obj))
     else:
         _baseclass = recarray
     attrdict = dict(_fieldmask=_fieldmask,
                     _hardmask=getattr(obj,'_hardmask',False),
                     _fill_value=getattr(obj,'_fill_value',None),
                     _sharedmask=getattr(obj,'_sharedmask',False),
                     _baseclass=_baseclass)
     self.__dict__.update(attrdict)
     # Finalize as a regular maskedarray .....
     # Update special attributes ...
     self._basedict = getattr(obj, '_basedict', getattr(obj,'__dict__',{}))
     self.__dict__.update(self._basedict)
     return

def mquantiles(data, prob=list([.25,.5,.75]), alphap=.4, betap=.4, axis=None):
    """Computes empirical quantiles for a *1xN* data array.
Samples quantile are defined by:
*Q(p) = (1-g).x[i] +g.x[i+1]*
where *x[j]* is the jth order statistic, 
with *i = (floor(n*p+m))*, *m=alpha+p*(1-alpha-beta)* and *g = n*p + m - i)*.

Typical values of (alpha,beta) are:
    
    - (0,1)    : *p(k) = k/n* : linear interpolation of cdf (R, type 4)
    - (.5,.5)  : *p(k) = (k+1/2.)/n* : piecewise linear function (R, type 5)
    - (0,0)    : *p(k) = k/(n+1)* : (R type 6)
    - (1,1)    : *p(k) = (k-1)/(n-1)*. In this case, p(k) = mode[F(x[k])].
      That's R default (R type 7)
    - (1/3,1/3): *p(k) = (k-1/3)/(n+1/3)*. Then p(k) ~ median[F(x[k])].
      The resulting quantile estimates are approximately median-unbiased
      regardless of the distribution of x. (R type 8)
    - (3/8,3/8): *p(k) = (k-3/8)/(n+1/4)*. Blom.
      The resulting quantile estimates are approximately unbiased
      if x is normally distributed (R type 9)
    - (.4,.4)  : approximately quantile unbiased (Cunnane)
    - (.35,.35): APL, used with PWM

:Parameters:
    x : Sequence
        Input data, as a sequence or array of dimension at most 2.
    prob : Sequence *[(0.25, 0.5, 0.75)]*
        List of quantiles to compute.
    alpha : Float (*[0.4]*)
        Plotting positions parameter.
    beta : Float (*[0.4]*)
        Plotting positions parameter.
    axis : Integer *[None]*
        Axis along which to compute quantiles. If *None*, uses the whole 
        (flattened/compressed) dataset.
    """
    def _quantiles1D(data,m,p):
        x = numpy.sort(data.compressed())
        n = len(x)
        if n == 0:
            return masked_array(numpy.empty(len(p), dtype=float_), mask=True)
        elif n == 1:
            return masked_array(numpy.resize(x, p.shape), mask=nomask)
        aleph = (n*p + m)
        k = numpy.floor(aleph.clip(1, n-1)).astype(int_)
        gamma = (aleph-k).clip(0,1)
        return (1.-gamma)*x[(k-1).tolist()] + gamma*x[k.tolist()]

    # Initialization & checks ---------
    data = masked_array(data, copy=False)
    p = narray(prob, copy=False, ndmin=1)
    m = alphap + p*(1.-alphap-betap)
    # Computes quantiles along axis (or globally)
    if (axis is None): 
        return _quantiles1D(data, m, p)
    else:
        assert data.ndim <= 2, "Array should be 2D at most !"
        return apply_along_axis(_quantiles1D, axis, data, m, p)

def cov(x, y=None, rowvar=True, bias=False, strict=False):
    """
    Estimate the covariance matrix.

    If x is a vector, return the variance.  For matrices, returns the covariance 
    matrix.

    If y is given, it is treated as an additional (set of) variable(s).

    Normalization is by (N-1) where N is the number of observations (unbiased 
    estimate).  If bias is True then normalization is by N.

    If rowvar is non-zero (default), then each row is a variable with observations 
    in the columns, otherwise each column is a variable  and the observations  are 
    in the rows.
    
    If strict is True, masked values are propagated: if a masked value appears in 
    a row or column, the whole row or column is considered masked.
    """
    X = narray(x, ndmin=2, subok=True, dtype=float)
    if X.shape[0] == 1:
        rowvar = True
    if rowvar:
        axis = 0
        tup = (slice(None),None)
    else:
        axis = 1
        tup = (None, slice(None))
    #
    if y is not None:
        y = narray(y, copy=False, ndmin=2, subok=True, dtype=float)
        X = concatenate((X,y),axis)
    #
    X -= X.mean(axis=1-axis)[tup]
    n = X.count(1-axis)
    #
    if bias:
        fact = n*1.0
    else:
        fact = n-1.0
    #
    if not rowvar:
        return (dot(X.T, X.conj(), strict=False) / fact).squeeze()
    else:
        return (dot(X, X.T.conj(), strict=False) / fact).squeeze()

 def write(self,file):
     iclose=False
     try:
         # If file is not open, open it as a fileobject
         fo=open(file,'wb')
         iclose=True
     except:
         # If already a fileobject, we will append to it
         fo=file
     narray(self.nboot).tofile(fo)
     narray(self.ndim).tofile(fo)
     self.Avg.tofile(fo)
     self.Std.tofile(fo)
     self.values.tofile(fo)
     
     # If file was a new file, close it
     if iclose:
         fo.close()

 def AddValue(self, value):
     """supply one more numerical value"""
     if(self.previous_index == self.size - 1
        and not(self.buffer_full)):
         self.buffer_full = True
     self.previous_index = (self.previous_index + 1) % self.size
     self.buf[self.previous_index * 2 + 1] = value
     self.vbo[self.previous_index * 2 + 1:
              self.previous_index * 2 + 2] = narray([value], 'f')
     self.maxValue = max(self.maxValue, value)

def plot_values(teamsdict, teamslist, vals):
    vals = narray(vals).T
    v = []
    for team in teamslist:
        if team in teamsdict:
            v.append(vals[teamsdict[team]])

    v = zip(*v)

    plt(v[0], v[1], 'rx')
    show()

 def _UpdateObject(self, _id, obj):
     """Updates an object in the VBO"""
     if self._vbo is None:
         return
     index = self._indices.get(_id, -1)
     # have to check since object could have been deleted since marked as
     # dirty
     if index < 0:
         return
     num_values = type(self).__num_values
     self._vbo[index * num_values : (index + 1) * num_values] = narray(self.__descToArray(obj), "f")

def hinton(W, length, name, 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))
		maxWeight = max(abs(W))*1.2

	P.fill(narray([0,width,width,0]),narray([0,0,height,height]),'gray')
	P.axis('off')
	P.axis('equal')
	correlationcounter = 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*w/maxWeight),'white')
				correlationcounter += w - 0.5
			elif w < 0:
				_blob(_x - 0.5, height - _y + 0.5, min(1,w/maxWeight*w/maxWeight),'black')
				correlationcounter += w + 0.5

		
	
	correlationindex = abs(correlationcounter/length*100 + 50)
	titlestring = 'My votes versus ' + name + ': ' + str(correlationindex) + '% agreement' 

	P.title(titlestring)
	if reenable:
		P.ion()
	P.show()

    def Import(self,data,bin=1,weights=[],randlist = []):
        if len(weights)>0 and bin!=1:
            print "Weights only currently supported for binsize=1"
            assert 1==0
        self.nconf=len(data)/bin
        if bin>1:
            tmpdata=narray([average(ien) for ien in split(narray(data[0:bin*(len(data)/bin)]),self.nconf)])
        else:
            tmpdata=narray(data)
#        self.Raw=append(self.Raw,tmpdata)
#        self.nconf=len(self.Raw)
#        rint=zeros( (self.nboot,self.nconf))
        # if randlist == []:
        myseed=startseed*self.nconf/self.nboot
        random.seed(myseed)
        #        locranint=random.randint
        locranint=random.random_integers
        # else:
        #     locranint = randlist
        if len(weights)>0:
            tmpweight=narray(weights)
            self.Avg=average(multiply(tmpdata,tmpweight))/sum(tmpweight)
        else:
            self.Avg=average(tmpdata)
        for iboot in range(self.nboot):
            rint=locranint(0,self.nconf-1,self.nconf)
#            tmp=0.0
#            for iconf in range(self.nconf):
#                rint=locranint(0,self.nconf-1)
#                tmp+=tmpdata[rint]/self.nconf
#            self.values[iboot]=tmp
            if len(weights)>0:
                tw2=tmpweight[rint]                
                td2=multiply(tmpdata[rint],tw2)
                self.values[iboot]=average(td2)/sum(tw2)
#average(multiply(tmpdata[rint],tmpweight[rint]))/sum(tmpweight[rint])
            else:
                self.values[iboot]=average(tmpdata[rint])
        return locranint

 def correlate(self):
     """ inst.correlate() -> None
     This must be overloaded to populate the '_correlations'
     member dictionary with valid correlation functions for 
     every included (i.e. tracked) baseline. This should be a
     numpy array.
     """
     with RLock():
         itime = self.server._integration_time
     self._stopevent.wait(itime)
     self._last_correlation = time()
     for baseline in self._include_baselines:
         self._correlations[baseline] = narray([0]*self._lags*2)

    def __getitem__(self, indx):
        """Returns all the fields sharing the same fieldname base.
The fieldname base is either `_data` or `_mask`."""
        _localdict = self.__dict__
        _fieldmask = _localdict['_fieldmask']
        _data = self._data
        # We want a field ........
        if isinstance(indx, basestring):
            obj = _data[indx].view(MaskedArray)
            obj._set_mask(_fieldmask[indx])
            # Force to nomask if the mask is empty
            if not obj._mask.any():
                obj._mask = nomask
            # Force to masked if the mask is True
            if not obj.ndim and obj._mask:
                return masked
            return obj
        # We want some elements ..
        # First, the data ........
        obj = narray(_data[indx], copy=False).view(mrecarray)
        obj._fieldmask = narray(_fieldmask[indx], copy=False).view(recarray)
        return obj

def dict_to_numpy(csv_data, columns_to_exclude = []):
    target = []
    data = []
    attributes = None

    for row in csv_data:
        target.append(row.pop('quality'))
        [row.pop(col) for col in columns_to_exclude]

        if attributes is None:
          attributes = row.keys() 

        data.append(row.values())

    target = narray(target)
    data = narray(data)

    return {
      'target': target,
      'data': data,
      'attributes': attributes
    }

 def __init__(self, size):
     self.size = size
     initialBuffer = []
     for x, y in enumerate(range(size)):
         initialBuffer.append(x)
         initialBuffer.append(y)
     self.buf = narray(initialBuffer, 'f')
     self.vbo = vbo.VBO(self.buf,
                        usage=GL_STREAM_DRAW,
                        target=GL_ARRAY_BUFFER,
                        size=4 * 2 * size)
     self.previous_index = -1
     self.buffer_full = False
     self.maxValue = 0

def MakePalette():
    "Generate a palette suitable for PIL.Image.ImageData.putpalette"
    return flatten((i,i,i) for i in xrange(256))
    
    # The below code mimicks the source from the original
    
    # If you want to use this, you need numpy. "easy_install numpy" or whatever 
    # for your distribution
    from numpy import array as narray
    
    # Some possible destination colours (sepia, or blueish, etc)
    #dst = narray([18,74,155])
    #dst = narray([155,74,18])
    #dst = narray([17,51,96])
    dst = narray([96,51,17])
    
    high = narray([255,255,255])
    firsthalf  = [high + (dst-high)/255*i for i in xrange(0,256,2)]
    secondhalf = [dst + (-dst)/255*i for i in xrange(0,256,2)]
    
    result = flatten(firsthalf+secondhalf)
    result = map(int, result)
    return result

def Plot4D( x, y, z, t, markerstyle = 20, markersize = 1 ):
    data = narray( [0.] * 4 )
    tree = TTree('DummyTree','DummyTree')
    tree.Branch('xyzt', data, 'x/D:y:z:t')
    
    for datai in zip(x,y,z,t):
        data[0], data[1], data[2], data[3] = datai
        #print data
        tree.Fill()
    tree.SetMarkerStyle( markerstyle )
    tree.SetMarkerSize( markersize )
    c = TCanvas()
    tree.Draw('x:y:z:t','','zcol')
    return c, tree

 def _AddObject(self, _id, obj):
     """Adds an object to the VBO"""
     if self._vbo is None:
         return
     index = -1
     if len(self._empty_indices) > 0:
         index = self._empty_indices.pop()
     elif self._max_index < self._data_size - 1:
         self._max_index += 1
         index = self._max_index
     if index > 0:
         values = self.__descToArray(obj)
         num_values = len(values)
         self._vbo[index * num_values : (index + 1) * num_values] = narray(values, "f")
         self._indices[_id] = index
     else:
         self._vbo = None