def _set_reach_dist(self, point_index, processed, X, nbrs):
        P = X[point_index:point_index + 1]
        # Assume that radius_neighbors is faster without distances
        # and we don't need all distances, nevertheless, this means
        # we may be doing some work twice.
        indices = nbrs.radius_neighbors(P, radius=self.max_eps,
                                        return_distance=False)[0]

        # Getting indices of neighbors that have not been processed
        unproc = np.compress((~np.take(processed, indices)).ravel(),
                             indices, axis=0)
        # Neighbors of current point are already processed.
        if not unproc.size:
            return

        # Only compute distances to unprocessed neighbors:
        if self.metric == 'precomputed':
            dists = X[point_index, unproc]
        else:
            dists = pairwise_distances(P, np.take(X, unproc, axis=0),
                                       self.metric, n_jobs=None).ravel()

        rdists = np.maximum(dists, self.core_distances_[point_index])
        improved = np.where(rdists < np.take(self.reachability_, unproc))
        self.reachability_[unproc[improved]] = rdists[improved]
        self.predecessor_[unproc[improved]] = point_index

 def gt_topk(dat, axis, ret_typ, k, is_ascend):
     if ret_typ == "indices":
         if is_ascend:
             indices = np.arange(k)
         else:
             indices = np.arange(-1, -k-1, -1)
         ret = np.take(dat.argsort(axis=axis), axis=axis, indices=indices, mode='wrap')
     elif ret_typ == "value":
         if is_ascend:
             indices = np.arange(k)
         else:
             indices = np.arange(-1, -k-1, -1)
         ret = np.take(np.sort(dat, axis=axis), axis=axis, indices=indices, mode='wrap')
     else:
         assert dat.shape == (5, 5, 5, 5)
         assert axis is None or axis ==1
         ret = np.zeros(dat.shape)
         if is_ascend:
             indices = np.arange(k)
         else:
             indices = np.arange(-1, -k-1, -1)
         gt_argsort = np.take(dat.argsort(axis=axis), axis=axis, indices=indices, mode='wrap')
         if axis is None:
             ret.ravel()[gt_argsort] = 1
         else:
             for i in range(5):
                 for j in range(5):
                     for k in range(5):
                         ret[i, gt_argsort[i, :, j, k], j, k] = 1
     return ret

def permute_2d(m, p):
    """Performs 2D permutation of matrix m according to p."""
    return m[p][:, p]
    # unused below
    m_t = np.transpose(m)
    r_t = np.take(m_t, p, axis=0)
    return np.take(np.transpose(r_t), p, axis=0)

 def convert_to_8_bit(self):
     """
     Convert 16-bit display data to 8-bit using a lookup table.
     """
     if self.intensity_scaling == 'autoscale':
         self.display_min = self.display_data_16.min()
         self.display_max = self.display_data_16.max()
         self._make_linear_lookup_table()
     elif self.intensity_scaling == 'median_filter_autoscale':
         filtered_image = ndimage.filters.median_filter(
             self.display_data_16, size=3, output=self.filtered_image)
         self.display_min = self.filtered_image.min()
         self.display_max = self.filtered_image.max()
         self._make_linear_lookup_table()
     if not hasattr(self, 'display_data_8'):
         self.display_data_8 = np.empty(
             self.buffer_shape[1:], dtype=np.uint8)
     np.take(self.lut, self.display_data_16, out=self.display_data_8)
     try:
         self.display_intensity_scaling_queue.get_nowait()
     except Queue.Empty:
         pass
     self.display_intensity_scaling_queue.put(
         (self.intensity_scaling, self.display_min, self.display_max))
     self.image = ArrayInterfaceImage(self.display_data_8, allow_copy=False)
     pyglet.gl.glTexParameteri( #Reset to no interpolation
             pyglet.gl.GL_TEXTURE_2D,
             pyglet.gl.GL_TEXTURE_MAG_FILTER,
             pyglet.gl.GL_NEAREST)
     if hasattr(self, 'window'):
         if not self.window.visible:
             self.window.set_visible(True)
     return None

def gather(mask):
    import glob

    flist = glob.glob(mask)
    for f in flist:
        print f
        dset = da.read_nc(f, ["csat", "lat", "cloudpts"])
        csat = dset["csat"].values
        lat = dset["lat"].values
        altitude = dset["csat"].altitude
        idx = dset["cloudpts"].values > 0
        del dset
        cpts = dict()
        nprof = dict()
        for l in lats:
            idx1 = np.where((lat >= lats[l][0]) & (lat < lats[l][1]))[0]
            idx2 = np.where((lat >= -lats[l][1]) & (lat < -lats[l][0]))[0]
            idx = np.concatenate([idx1, idx2])
            if l in cpts:
                nprof[l] = nprof[l] + idx.shape[0]
                cpts[l] = cpts[l] + np.take(csat, idx, axis=0).sum(axis=0)
            else:
                nprof[l] = np.sum(idx)
                cpts[l] = np.take(csat, idx, axis=0).sum(axis=0)
    cprofl = dict()
    for l in lats:
        cprofl[l] = 100.0 * cpts[l] / nprof[l]

    return cprofl, altitude

def reconstruct(efds, T, K):
    
    T=np.ceil(T)
    
    N=len(efds)
    
    reconstructed = np.zeros((T,2))
    
    n = np.arange(start=1,stop=N,step=1)
    t = np.arange(T)
    
    n_grid, t_grid = np.meshgrid( n, t )
    
    a_n_grid = np.take(efds[:,0], n_grid)
    b_n_grid = np.take(efds[:,1], n_grid)
    c_n_grid = np.take(efds[:,2], n_grid)
    d_n_grid = np.take(efds[:,3], n_grid)

    arg_grid = n_grid * t_grid / T
    
    cos_term = np.cos( 2 * np.pi * arg_grid )
    sin_term = np.sin( 2 * np.pi * arg_grid )
    
    reconstructed[:,0] = efds[0,0] + np.sum(a_n_grid * cos_term + b_n_grid * sin_term, axis=1)
    reconstructed[:,1] = efds[0,0] + np.sum(c_n_grid * cos_term + d_n_grid * sin_term, axis=1)

    return reconstructed

def is_quadrant_red(quadrant):
     indices = xrange(0,len(quadrant),3)
     red_quadrant = np.take(quadrant, indices)
     red_sum = np.sum(red_quadrant)
     red_avg = np.sum(red_quadrant) / len(red_quadrant)

     logging.debug("red avg: %s" % (red_avg))

     indices = xrange(1,len(quadrant),3)
     green_quadrant = np.take(quadrant, indices)
     green_sum = np.sum(green_quadrant)
     green_avg = np.sum(green_quadrant) / len(green_quadrant)

     logging.debug("green avg: %s" % (green_avg))

     indices = xrange(2,len(quadrant),3)
     blue_quadrant = np.take(quadrant, indices)
     blue_sum = np.sum(blue_quadrant)
     blue_avg = np.sum(blue_quadrant) / len(blue_quadrant)

     logging.debug("blue avg: %s" % (blue_avg))

     is_red = red_avg / (0.5 * (green_avg + blue_avg))
     logging.debug("redcalc: %s" % (is_red))

     if is_red > 2:
         return 1
     else:
         return 0

    def interpolate(self, points):
        if self.tri == None:
            xc = self.x_coords.flatten()
            yc = self.y_coords.flatten()
            self.no_nan_values = self.values.flatten()

            if np.isnan(xc).any() and np.isnan(yc).any():
                xc = xc[~np.isnan(xc)]
                yc = yc[~np.isnan(yc)]
                self.no_nan_values = self.no_nan_values[~np.isnan(self.no_nan_values)]

            # Default: Qbb Qc Qz 
            self.tri = qhull.Delaunay(np.column_stack((xc, yc)), qhull_options='QbB')

        simplices = self.tri.find_simplex(points)

        indices = np.take(self.tri.simplices, simplices, axis=0)
        transforms = np.take(self.tri.transform, simplices, axis=0)

        delta = points - transforms[:,2]
        bary = np.einsum('njk,nk->nj', transforms[:,:2,:], delta)

        temp = np.hstack((bary, 1-bary.sum(axis=1, keepdims=True)))

        values = np.einsum('nj,nj->n', np.take(self.no_nan_values, indices), temp)

        #print values[np.any(temp<0, axis=1)]

        # This should put a NaN for points outside of any simplices
        # but is for some reason sometimes also true inside a simplex
        #values[np.any(temp < 0.0, axis=1)] = np.nan

        return values

    def __init__(self, x, y, ival=0., sorted=False, side='left'):

        if side.lower() not in ['right', 'left']:
            msg = "side can take the values 'right' or 'left'"
            raise ValueError(msg)
        self.side = side

        _x = np.asarray(x)
        _y = np.asarray(y)

        if _x.shape != _y.shape:
            msg = "x and y do not have the same shape"
            raise ValueError(msg)
        if len(_x.shape) != 1:
            msg = 'x and y must be 1-dimensional'
            raise ValueError(msg)

        self.x = np.r_[-np.inf, _x]
        self.y = np.r_[ival, _y]

        if not sorted:
            asort = np.argsort(self.x)
            self.x = np.take(self.x, asort, 0)
            self.y = np.take(self.y, asort, 0)
        self.n = self.x.shape[0]

def _set_reach_dist(core_distances_, reachability_, predecessor_,
                    point_index, processed, X, nbrs, metric, metric_params,
                    p, max_eps):
    P = X[point_index:point_index + 1]
    # Assume that radius_neighbors is faster without distances
    # and we don't need all distances, nevertheless, this means
    # we may be doing some work twice.
    indices = nbrs.radius_neighbors(P, radius=max_eps,
                                    return_distance=False)[0]

    # Getting indices of neighbors that have not been processed
    unproc = np.compress(~np.take(processed, indices), indices)
    # Neighbors of current point are already processed.
    if not unproc.size:
        return

    # Only compute distances to unprocessed neighbors:
    if metric == 'precomputed':
        dists = X[point_index, unproc]
    else:
        _params = dict() if metric_params is None else metric_params.copy()
        if metric == 'minkowski' and 'p' not in _params:
            # the same logic as neighbors, p is ignored if explicitly set
            # in the dict params
            _params['p'] = p
        dists = pairwise_distances(P, np.take(X, unproc, axis=0),
                                   metric, n_jobs=None,
                                   **_params).ravel()

    rdists = np.maximum(dists, core_distances_[point_index])
    improved = np.where(rdists < np.take(reachability_, unproc))
    reachability_[unproc[improved]] = rdists[improved]
    predecessor_[unproc[improved]] = point_index

    def _map(self, X):
        """ Maps from a scalar or an array to an RGBA value or array.

        The *X* parameter is either a scalar or an array (of any dimension).
        If it is scalar, the function returns a tuple of RGBA values; otherwise
        it returns an array with the new shape = oldshape+(4,).  Any values
        that are outside the 0,1 interval are clipped to that interval before
        generating RGB values.
        """

        if type(X) in [IntType, FloatType]:
            vtype = 'scalar'
            xa = array([X])
        else:
            vtype = 'array'
            xa = asarray(X)

        # assume the data is properly normalized
        #xa = where(xa>1.,1.,xa)
        #xa = where(xa<0.,0.,xa)


        nanmask = isnan(xa)
        xa = where(nanmask, 0, (xa * (self.steps-1)).astype(int))
        rgba = zeros(xa.shape+(4,), float)
        rgba[...,0] = where(nanmask, 0, take(self._red_lut, xa))
        rgba[...,1] = where(nanmask, 0, take(self._green_lut, xa))
        rgba[...,2] = where(nanmask, 0, take(self._blue_lut, xa))
        rgba[...,3] = where(nanmask, 0, take(self._alpha_lut, xa))
        if vtype == 'scalar':
            rgba = tuple(rgba[0,:])

        return rgba

    def test_np_ufuncs(self):
        z = self.create_array(shape=(100, 100), chunks=(10, 10))
        a = np.arange(10000).reshape(100, 100)
        z[:] = a

        eq(np.sum(a), np.sum(z))
        assert_array_equal(np.sum(a, axis=0), np.sum(z, axis=0))
        eq(np.mean(a), np.mean(z))
        assert_array_equal(np.mean(a, axis=1), np.mean(z, axis=1))
        condition = np.random.randint(0, 2, size=100, dtype=bool)
        assert_array_equal(np.compress(condition, a, axis=0),
                           np.compress(condition, z, axis=0))
        indices = np.random.choice(100, size=50, replace=True)
        assert_array_equal(np.take(a, indices, axis=1),
                           np.take(z, indices, axis=1))

        # use zarr array as indices or condition
        zc = self.create_array(shape=condition.shape, dtype=condition.dtype,
                               chunks=10, filters=None)
        zc[:] = condition
        assert_array_equal(np.compress(condition, a, axis=0),
                           np.compress(zc, a, axis=0))
        zi = self.create_array(shape=indices.shape, dtype=indices.dtype,
                               chunks=10, filters=None)
        zi[:] = indices
        # this triggers __array__() call with dtype argument
        assert_array_equal(np.take(a, indices, axis=1),
                           np.take(a, zi, axis=1))

    def _set_reach_dist(self, point_index, X, nbrs):
        P = np.array(X[point_index]).reshape(1, -1)
        indices = nbrs.radius_neighbors(P, radius=self.max_bound,
                                        return_distance=False)[0]

        # Getting indices of neighbors that have not been processed
        unproc = np.compress((~np.take(self._processed, indices)).ravel(),
                             indices, axis=0)
        # Keep n_jobs = 1 in the following lines...please
        if len(unproc) > 0:
            dists = pairwise_distances(P, np.take(X, unproc, axis=0),
                                       self.metric, n_jobs=1).ravel()

            rdists = np.maximum(dists, self.core_distances_[point_index])
            new_reach = np.minimum(np.take(self.reachability_, unproc), rdists)
            self.reachability_[unproc] = new_reach

        # Checks to see if everything is already processed;
        # if so, return control to main loop
        if unproc.size > 0:
            # Define return order based on reachability distance
            return(unproc[quick_scan(np.take(self.reachability_, unproc),
                                     dists)])
        else:
            return point_index

    def get_left_channels(self, energy, nchan=1):
        self.initialize()
        g_s_ii = self.greenfunction.retarded(energy)
        lambda_l_ii = self.selfenergies[0].get_lambda(energy)
        lambda_r_ii = self.selfenergies[1].get_lambda(energy)

        if self.greenfunction.S is not None:
            s_mm = self.greenfunction.S
            s_s_i, s_s_ii = linalg.eig(s_mm)
            s_s_i = np.abs(s_s_i)
            s_s_sqrt_i = np.sqrt(s_s_i)  # sqrt of eigenvalues
            s_s_sqrt_ii = np.dot(s_s_ii * s_s_sqrt_i, dagger(s_s_ii))
            s_s_isqrt_ii = np.dot(s_s_ii / s_s_sqrt_i, dagger(s_s_ii))

        lambdab_r_ii = np.dot(np.dot(s_s_isqrt_ii, lambda_r_ii), s_s_isqrt_ii)
        a_l_ii = np.dot(np.dot(g_s_ii, lambda_l_ii), dagger(g_s_ii))
        ab_l_ii = np.dot(np.dot(s_s_sqrt_ii, a_l_ii), s_s_sqrt_ii)
        lambda_i, u_ii = linalg.eig(ab_l_ii)
        ut_ii = np.sqrt(lambda_i / (2.0 * np.pi)) * u_ii
        m_ii = 2 * np.pi * np.dot(np.dot(dagger(ut_ii), lambdab_r_ii), ut_ii)
        T_i, c_in = linalg.eig(m_ii)
        T_i = np.abs(T_i)

        channels = np.argsort(-T_i)[:nchan]
        c_in = np.take(c_in, channels, axis=1)
        T_n = np.take(T_i, channels)
        v_in = np.dot(np.dot(s_s_isqrt_ii, ut_ii), c_in)

        return T_n, v_in

 def get_MW(self, F, mode='F^-1'):
     if type(F) is dict: # recursive case for many F's at once
         M,W = {}, {}
         for key in F: M[key],W[key] = self.get_MW(F[key], mode=mode)
         return M,W
     modes = ['F^-1', 'F^-1/2', 'I', 'L^-1']; assert(mode in modes)
     if mode == 'F^-1':
         M = np.linalg.pinv(F, rcond=1e-12)
         #U,S,V = np.linalg.svd(F)
         #M = np.einsum('ij,j,jk', V.T, 1./S, U.T)
     elif mode == 'F^-1/2':
         U,S,V = np.linalg.svd(F)
         M = np.einsum('ij,j,jk', V.T, 1./np.sqrt(S), U.T)
     elif mode == 'I':
         M = np.identity(F.shape[0], dtype=F.dtype)
     else:
         #Cholesky decomposition to get M
         order = np.array([10,11,9,12,8,20,0,13,7,14,6,15,5,16,4,17,3,18,2,19,1]) # XXX needs generalizing
         iorder = np.argsort(order)
         F_o = np.take(np.take(F,order, axis=0), order, axis=1)
         L_o = np.linalg.cholesky(F_o)
         U,S,V = np.linalg.svd(L_o.conj())
         M_o = np.dot(np.transpose(V), np.dot(np.diag(1./S), np.transpose(U)))
         M = np.take(np.take(M_o,iorder, axis=0), iorder, axis=1)
     W = np.dot(M, F)
     norm  = W.sum(axis=-1); norm.shape += (1,)
     M /= norm; W = np.dot(M, F)
     return M,W

def plot_eig(data,nchan):
    days=data.keys()
    for k in days:
        eig_order=[]
        eigs = []
        eigs_cav = []
        for bl in data[k]:
           c_mat=cov(data[k][bl])
           cav = get_cav(c_mat,nchan,scaling=opts.auto)
           U,S,V= n.linalg.svd(c_mat.conj())
           U_cav,S_cav,V_cav = n.linalg.svd(cav.conj())
           eig_order.append(S[0])
           eigs.append( n.fft.fftshift(n.fft.fft(V.T.conj(),axis=0)))
           eigs_cav.append( n.fft.fftshift(n.fft.fft(V_cav.T.conj(),axis=0)))


        order=n.argsort(eig_order)

        eig_order=n.take(eig_order,order)
        eigs=n.take(eigs,order,axis=0)
        eigs_cav=n.take(eigs_cav,order,axis=0)
        embed()
        fig=p.figure(1)
        for cnt,eig in enumerate(eigs):
            p.plot(eig[0] + cnt*5)
        p.title('Eigenvectors for day {0}'.format(k))
        p.show()
        p.savefig('eigenvectors_{0}.png'.format(k))
        p.clf()
        for cnt,eig in enumerate(eigs_cav):
            p.plot(eig[0] + cnt*5)
        p.title('Eigenvectors of Cav for day {0}'.format(k))
        p.savefig('eigenvectors_cav_{0}.png'.format(k))
        p.clf()
        p.close()

def continuous_components(delta_X, delta_Y, delta_t, t, T, K):
       
    p = np.arange(K)
    
    delta_xp = np.take(delta_X, p)
    delta_yp = np.take(delta_Y, p)
    delta_tp = np.take(delta_t, p)
    tp = np.take(t, p)
    tp = np.hstack( ( np.array([0]) , tp ) )
    
    first_term_xi = np.cumsum(delta_X[0:K-1])
    second_term_xi = (delta_X[1:K]/delta_t[1:K]) * np.cumsum(delta_t[0:K-1])
    xi = np.hstack( ( np.array([0]), first_term_xi - second_term_xi ) )
    
    first_term_delta = np.cumsum(delta_Y[0:K-1])
    second_term_delta = (delta_Y[1:K]/delta_t[1:K]) * np.cumsum(delta_t[0:K-1])
    delta = np.hstack( ( np.array([0]), first_term_delta - second_term_delta ) )
        
    A0 = (1/T)*np.sum( (delta_xp/(2*delta_tp) * (np.square(tp[1:K+1]) - np.square(tp[0:K]))) + \
                       xi * (tp[1:K+1] - tp[0:K]))

    C0 = (1/T)*np.sum( (delta_yp/(2*delta_tp) * (np.square(tp[1:K+1]) - np.square(tp[0:K]))) + \
                       delta * (tp[1:K+1] - tp[0:K]))
    
    return A0, C0

	def resample(self):
		"resample() randomly draws a set of points equal in size to the original set from the cached data for bootstrapping"
		assert hasattr(self, "saved_xarray"), "resampling not set up yet.  Call setup_resampling() first."
		ranlist=Numeric.floor(self.get_random_list(self.pointcount)*self.pointcount).astype(numeric_int)
		self.xarray=Numeric.take(self.saved_xarray, ranlist, -1) #take columns since vectors lie this way
		self.yarray=Numeric.take(self.saved_yarray, ranlist)
		self.firstpass=1

 def test_take_output(self, level=rlevel):
     """Ensure that 'take' honours output parameter."""
     x = np.arange(12).reshape((3,4))
     a = np.take(x,[0,2],axis=1)
     b = np.zeros_like(a)
     np.take(x,[0,2],axis=1,out=b)
     assert_array_equal(a,b)

def shoelace(vertices):
    """        
    Calculate twice the area of polygon using Shoelace formula.
    
    Polygon is defined by vertices.
    
    Parameters
    ----------
    vertices : array_like
        Vertex coordinates in a 2-D space.
        Coordinates must be placed along the last axis. And data points are
        along the first axis.
    
    Returns
    -------
    area : float
        You can deduce the order of input vertices from the sign:
        area is positive if vertices are in counter-clockwise order.
        area is negative if vertices are in clockwise order.
        area is zero if all points are colinear.
    
    Notes
    -----
    This function can be also used to judge if all points in a data set are
    collinear. Collinear points as input for initializing Polygon instance
    will raise a QhullError.
    
    Examples
    --------
    Vertices of a square:

    Clockwise:    
    
    >>> from tadlib.calfea.polygon import shoelace
    >>> sq = [(0,0), (0,1), (1,1), (1,0)]
    >>> shoelace(sq)
    -2.0
    
    Counter-clockwise:
    
    >>> sq = [(0,0), (1,0), (1,1), (0,1)]
    >>> shoelace(sq)
    2.0
        
    """
    
    vertices = np.asfarray(vertices)
    # Rule for stacking multiple comma separated arrays
    rule = '0,' + str(len(vertices.shape))
    # Slip the array along the first axis
    slip_v = np.r_[rule, vertices[-1], vertices[:-1]]
    # Extract coordinates
    x = np.take(vertices, [0], axis=-1).reshape(vertices.shape[:-1])
    y = np.take(vertices, [1], axis=-1).reshape(vertices.shape[:-1])
    slip_x = np.take(slip_v, [0], axis=-1).reshape(vertices.shape[:-1])
    slip_y = np.take(slip_v, [1], axis=-1).reshape(vertices.shape[:-1])
    # Sholelace Foluma
    area = np.sum(y * slip_x - x * slip_y, axis=0)
    
    return area