How to use the healpy.get_all_neighbours function in healpy

npix = len(m)
    nside = hp.npix2nside(npix)
    min_area = 0.4 * hp.nside2pixarea(nside)

    # Compute faces, vertices, and neighbors.
    # vertices is an N X 3 array of the distinct vertices of the HEALPix faces.
    # faces is an npix X 4 array mapping HEALPix faces to their vertices.
    # neighbors is an npix X 4 array mapping faces to their nearest neighbors.
    faces = np.ascontiguousarray(
        np.rollaxis(hp.boundaries(nside, np.arange(npix), nest=nest), 2, 1))
    dtype = faces.dtype
    faces = faces.view(np.dtype((np.void, dtype.itemsize * 3)))
    vertices, faces = np.unique(faces.ravel(), return_inverse=True)
    faces = faces.reshape(-1, 4)
    vertices = vertices.view(dtype).reshape(-1, 3)
    neighbors = hp.get_all_neighbours(nside, np.arange(npix), nest=nest)[::2].T

    # Loop over the requested contours.
    paths = []
    for level in levels:

        # Find credible region
        indicator = (m >= level)

        # Construct a graph of the eges of the contour.
        graph = nx.Graph()
        face_pairs = set()
        for ipix1, ipix2 in enumerate(neighbors):
            for ipix2 in ipix2:
                # Determine if we have already considered this pair of faces.
                new_face_pair = frozenset((ipix1, ipix2))
                if new_face_pair in face_pairs:

def flood_fill(nside, ipix, m, nest=False):
    """Stack-based flood fill algorithm in HEALPix coordinates.
    Based on .
    """
    # Initialize stack with starting pixel index.
    stack = [ipix]
    while stack:
        # Pop last pixel off of the stack.
        ipix = stack.pop()
        # Is this pixel in need of filling?
        if m[ipix]:
            # Fill in this pixel.
            m[ipix] = False
            # Find the pixels neighbors.
            neighbors = hp.get_all_neighbours(nside, ipix, nest=nest)
            # All pixels have up to 8 neighbors. If a pixel has less than 8
            # neighbors, then some entries of the array are set to -1. We
            # have to skip those.
            neighbors = neighbors[neighbors != -1]
            # Push neighboring pixels onto the stack.
            stack.extend(neighbors)

-------
        List of observations.
        """
        obs_array = np.concatenate(observation_list)
        hpids = _raDec2Hpid(self.nside, obs_array['RA'], obs_array['dec'])
        new_expts = np.zeros(obs_array.size, dtype=float)
        for filtername in np.unique(obs_array['filter']):
            in_filt = np.where(obs_array['filter'] == filtername)
            delta_m5 = self.target_m5[filtername] - conditions.M5Depth[filtername][hpids[in_filt]]
            # We can get NaNs because dithering pushes the center of the pointing into masked regions.
            nan_indices = np.argwhere(np.isnan(delta_m5)).ravel()
            for indx in nan_indices:
                bad_hp = hpids[in_filt][indx]
                # Note this might fail if we run at higher resolution, then we'd need to look farther for
                # pixels to interpolate.
                near_pix = hp.get_all_neighbours(conditions.nside, bad_hp)
                vals = conditions.M5Depth[filtername][near_pix]
                if True in np.isfinite(vals):
                    estimate_m5 = np.mean(vals[np.isfinite(vals)])
                    delta_m5[indx] = self.target_m5[filtername] - estimate_m5
                else:
                    raise ValueError('Failed to find a nearby unmasked sky value.')

            new_expts[in_filt] = conditions.exptime * 10**(delta_m5/1.25)
        new_expts = np.clip(new_expts, self.min_exp, self.max_exp)
        # I'm not sure what level of precision we can expect, so let's just limit to seconds
        new_expts = np.round(new_expts)

        for i, observation in enumerate(observation_list):
            observation['exptime'] = new_expts[i]

        return observation_list

def _get_neighbors(self, idx):
        import healpy as hp

        nside = self._get_nside(idx)
        idx_nb = (hp.get_all_neighbours(nside, idx[0], nest=self.nest),)
        idx_nb += tuple([t[None, ...] * np.ones_like(idx_nb[0]) for t in idx[1:]])

        return idx_nb

def flood_fill(nside, ipix, m, nest=False):
    """Stack-based flood fill algorithm in HEALPix coordinates.
    Based on .
    """
    # Initialize stack with starting pixel index.
    stack = [ipix]
    while stack:
        # Pop last pixel off of the stack.
        ipix = stack.pop()
        # Is this pixel in need of filling?
        if m[ipix]:
            # Fill in this pixel.
            m[ipix] = False
            # Find the pixels neighbors.
            neighbors = hp.get_all_neighbours(nside, ipix, nest=nest)
            # All pixels have up to 8 neighbors. If a pixel has less than 8
            # neighbors, then some entries of the array are set to -1. We
            # have to skip those.
            neighbors = neighbors[neighbors != -1]
            # Push neighboring pixels onto the stack.
            stack.extend(neighbors)

def analysis(targ_ra, targ_dec, mod):
    """Analyze a candidate"""

    #hdisc = healpix.ang2disc(nside, targ_ra, targ_dec, 1, inclusive=True)
    #files = []
    #for i in hdisc:
    #    if os.path.exists(datadir+'/cat_hpx_'+str(i).zfill(5)+'.fits'):
    #        files.append(datadir+'/cat_hpx_'+str(i).zfill(5)+'.fits')
    #
    ##data = utils.load_infiles(files, columns=['RA', 'DEC', 'WAVG_MAG_PSF_G', 'WAVG_MAG_PSF_R', 'EXTINCTION_G', 'EXTINCTION_R', 'WAVG_SPREAD_MODEL_R', 'SPREADERR_MODEL_R', 'WAVG_MAGERR_PSF_G', 'WAVG_MAGERR_PSF_R', 'MAGERR_PSF_G', 'MAGERR_PSF_R'])
    #data = utils.load_infiles(files, columns=COLUMNS)

    pix_nside_select = ugali.utils.healpix.angToPix(nside, targ_ra, targ_dec)
    ra_select, dec_select = ugali.utils.healpix.pixToAng(nside, pix_nside_select)
    pix_nside_neighbors = np.concatenate([[pix_nside_select], healpy.get_all_neighbours(nside, pix_nside_select)])
    data_array = []
    for pix_nside in pix_nside_neighbors:
        #infile = '%s/cat_hpx_%05i.fits'%(datadir, pix_nside)
        infile = '%s/y3a2_ngmix_cm_%05i.fits'%(datadir, pix_nside)
        #infile = '%s/*_%05i.fits'%(datadir, pix_nside) # TODO: get to work
        if not os.path.exists(infile):
            continue
        reader = pyfits.open(infile)
        data_array.append(reader[1].data)
        reader.close()
    print('Assembling data...')
    data = np.concatenate(data_array)
    print('Found {} objects...').format(len(data))
    print('Loading data...')

    ## De-redden magnitudes