docs: def randomPositions(input, nside_pix, n=1):
      """
      Generate n random positions within a full HEALPix mask of booleans, or a set of (lon, lat) coordinates.
  
      Parameters:
      -----------
      input :     (1) full HEALPix mask of booleans, or (2) a set of (lon, lat) coordinates for catalog objects that define the occupied pixels.
      nside_pix : nside_pix is meant to be at coarser resolution than the input mask or catalog object positions
      so that gaps from star holes, bleed trails, cosmic rays, etc. are filled in. 
  
      Returns:
      --------
      lon,lat,area : Return the longitude and latitude of the random positions (deg) and the total area (deg^2).
  
      """
      input = np.array(input)
      if len(input.shape) == 1:
          if hp.npix2nside(len(input)) < nside_pix:
              logger.warning('Expected coarser resolution nside_pix in skymap.randomPositions')
          subpix = np.nonzero(input)[0] # All the valid pixels in the mask at the NSIDE for the input mask
          lon, lat = pix2ang(hp.npix2nside(len(input)), subpix)
      elif len(input.shape) == 2:
          lon, lat = input[0], input[1] # All catalog object positions
      else:
          logger.warning('Unexpected input dimensions for skymap.randomPositions')
      pix = surveyPixel(lon, lat, nside_pix)
  
      # Area with which the random points are thrown
      area = len(pix) * hp.nside2pixarea(nside_pix, degrees=True)
  
      # Create mask at the coarser resolution
      mask = np.tile(False, hp.nside2npix(nside_pix))
      mask[pix] = True
  
      # Estimate the number of points that need to be thrown based off
      # coverage fraction of the HEALPix mask
      coverage_fraction = float(np.sum(mask)) / len(mask) 
      n_throw = int(n / coverage_fraction)
          
      lon, lat = [], []
      count = 0
      while len(lon) < n:
          lon_throw = np.random.uniform(0., 360., n_throw)
          lat_throw = np.degrees(np.arcsin(np.random.uniform(-1., 1., n_throw)))
  
          pix_throw = ugali.utils.healpix.angToPix(nside_pix, lon_throw, lat_throw)
          cut = mask[pix_throw].astype(bool)
  
          lon = np.append(lon, lon_throw[cut])
          lat = np.append(lat, lat_throw[cut])
  
          count += 1
          if count > 10:
              raise RuntimeError('Too many loops...')
  
      return lon[0:n], lat[0:n], area