docs: def query(self, coords, mode='random_sample', return_flags=False, pct=None):
          """
          Returns reddening at the requested coordinates. There are several
          different query modes, which handle the probabilistic nature of the map
          differently.
  
          Args:
              coords (:obj:`astropy.coordinates.SkyCoord`): The coordinates to query.
              mode (Optional[:obj:`str`]): Seven different query modes are available:
                  'random_sample', 'random_sample_per_pix' 'samples', 'median',
                  'mean', 'best' and 'percentile'. The :obj:`mode` determines how the
                  output will reflect the probabilistic nature of the Bayestar
                  dust maps.
              return_flags (Optional[:obj:`bool`]): If :obj:`True`, then QA flags will be
                  returned in a second numpy structured array. That is, the query
                  will return :obj:`ret`, :obj:'flags`, where :obj:`ret` is the normal return
                  value, containing reddening. Defaults to :obj:`False`.
              pct (Optional[:obj:`float` or list/array of :obj:`float`]): If the mode is
                  :obj:`percentile`, then :obj:`pct` specifies which percentile(s) is
                  (are) returned.
  
          Returns:
              Reddening at the specified coordinates, in magnitudes of reddening.
  
              The conversion to E(B-V) (or other reddening units) depends on
              whether :obj:`version='bayestar2017'` (the default) or
              :obj:`'bayestar2015'` was selected when the :obj:`BayestarQuery` object
              was created. To convert Bayestar2017 to Pan-STARRS 1 extinctions,
              multiply by the coefficients given in Table 1 of Green et al.
              (2018). Conversion to extinction in non-PS1 passbands depends on the
              choice of extinction law. To convert Bayestar2015 to extinction in
              various passbands, multiply by the coefficients in Table 6 of
              Schlafly & Finkbeiner (2011). See Green et al. (2015, 2018) for more
              detailed discussion of how to convert the Bayestar dust maps into
              reddenings or extinctions in different passbands.
  
              The shape of the output depends on the :obj:`mode`, and on whether
              :obj:`coords` contains distances.
  
              If :obj:`coords` does not specify distance(s), then the shape of the
              output begins with :obj:`coords.shape`. If :obj:`coords` does specify
              distance(s), then the shape of the output begins with
              :obj:`coords.shape + ([number of distance bins],)`.
  
              If :obj:`mode` is :obj:`'random_sample'`, then at each
              coordinate/distance, a random sample of reddening is given.
  
              If :obj:`mode` is :obj:`'random_sample_per_pix'`, then the sample chosen
              for each angular pixel of the map will be consistent. For example,
              if two query coordinates lie in the same map pixel, then the same
              random sample will be chosen from the map for both query
              coordinates.
  
              If :obj:`mode` is :obj:`'median'`, then at each coordinate/distance, the
              median reddening is returned.
  
              If :obj:`mode` is :obj:`'mean'`, then at each coordinate/distance, the
              mean reddening is returned.
  
              If :obj:`mode` is :obj:`'best'`, then at each coordinate/distance, the
              maximum posterior density reddening is returned (the "best fit").
  
              If :obj:`mode` is :obj:`'percentile'`, then an additional keyword
              argument, :obj:`pct`, must be specified. At each coordinate/distance,
              the requested percentiles (in :obj:`pct`) will be returned. If :obj:`pct`
              is a list/array, then the last axis of the output will correspond to
              different percentiles.
  
              Finally, if :obj:`mode` is :obj:`'samples'`, then at each
              coordinate/distance, all samples are returned. The last axis of the
              output will correspond to different samples.
  
              If :obj:`return_flags` is :obj:`True`, then in addition to reddening, a
              structured array containing QA flags will be returned. If the input
              coordinates include distances, the QA flags will be :obj:`"converged"`
              (whether or not the line-of-sight fit converged in a given pixel)
              and :obj:`"reliable_dist"` (whether or not the requested distance is
              within the range considered reliable, based on the inferred
              stellar distances). If the input coordinates do not include
              distances, then instead of :obj:`"reliable_dist"`, the flags will
              include :obj:`"min_reliable_distmod"` and :obj:`"max_reliable_distmod"`,
              the minimum and maximum reliable distance moduli in the given pixel.
          """
  
          # Check that the query mode is supported
          self._raise_on_mode(mode)
  
          # Validate percentile specification
          pct, scalar_pct = self._interpret_percentile(mode, pct)
  
          # Get number of coordinates requested
          n_coords_ret = coords.shape[0]
  
          # Determine if distance has been requested
          has_dist = hasattr(coords.distance, 'kpc')
          d = coords.distance.kpc if has_dist else None
  
          # Extract the correct angular pixel(s)
          # t0 = time.time()
          pix_idx = self._find_data_idx(coords.l.deg, coords.b.deg)
          in_bounds_idx = (pix_idx != -1)
  
          # t1 = time.time()
  
          # Extract the correct samples
          if mode == 'random_sample':
              # A different sample in each queried coordinate
              samp_idx = np.random.randint(0, self._n_samples, pix_idx.size)
              n_samp_ret = 1
          elif mode == 'random_sample_per_pix':
              # Choose same sample in all coordinates that fall in same angular
              # HEALPix pixel
              samp_idx = np.random.randint(0, self._n_samples, self._n_pix)[pix_idx]
              n_samp_ret = 1
          elif mode == 'best':
              samp_idx = slice(None)
              n_samp_ret = 1
          else:
              # Return all samples in each queried coordinate
              samp_idx = slice(None)
              n_samp_ret = self._n_samples
  
          # t2 = time.time()
  
          if mode == 'best':
              val = self._best_fit
          else:
              val = self._samples
  
          # Create empty array to store flags
          if return_flags:
              if has_dist:
                  # If distances are provided in query, return only covergence and
                  # whether or not this distance is reliable
                  dtype = [('converged', 'bool'),
                           ('reliable_dist', 'bool')]
                  # shape = (n_coords_ret)
              else:
                  # Return convergence and reliable distance ranges
                  dtype = [('converged', 'bool'),
                           ('min_reliable_distmod', 'f4'),
                           ('max_reliable_distmod', 'f4')]
              flags = np.empty(n_coords_ret, dtype=dtype)
          # samples = self._samples[pix_idx, samp_idx]
          # samples[pix_idx == -1] = np.nan
  
          # t3 = time.time()
  
          # Extract the correct distance bin (possibly using linear interpolation)
          if has_dist: # Distance has been provided
              # Determine ceiling bin index for each coordinate
              dm = 5. * (np.log10(d) + 2.)
              bin_idx_ceil = np.searchsorted(self._DM_bin_edges, dm)
  
              # Create NaN-filled return arrays
              if isinstance(samp_idx, slice):
                  ret = np.full((n_coords_ret, n_samp_ret), np.nan, dtype='f4')
              else:
                  ret = np.full((n_coords_ret,), np.nan, dtype='f4')
  
              # d < d(nearest distance slice)
              idx_near = (bin_idx_ceil == 0) & in_bounds_idx
              if np.any(idx_near):
                  a = 10.**(0.2 * (dm[idx_near] - self._DM_bin_edges[0]))
                  if isinstance(samp_idx, slice):
                      ret[idx_near] = (
                          a[:,None]
                          * val[pix_idx[idx_near], samp_idx, 0])
                  else:
                      # print('idx_near: {} true'.format(np.sum(idx_near)))
                      # print('ret[idx_near].shape = {}'.format(ret[idx_near].shape))
                      # print('val.shape = {}'.format(val.shape))
                      # print('pix_idx[idx_near].shape = {}'.format(pix_idx[idx_near].shape))
  
                      ret[idx_near] = (
                          a * val[pix_idx[idx_near], samp_idx[idx_near], 0])
  
              # d > d(farthest distance slice)
              idx_far = (bin_idx_ceil == self._n_distances) & in_bounds_idx
              if np.any(idx_far):
                  # print('idx_far: {} true'.format(np.sum(idx_far)))
                  # print('pix_idx[idx_far].shape = {}'.format(pix_idx[idx_far].shape))
                  # print('ret[idx_far].shape = {}'.format(ret[idx_far].shape))
                  # print('val.shape = {}'.format(val.shape))
                  if isinstance(samp_idx, slice):
                      ret[idx_far] = val[pix_idx[idx_far], samp_idx, -1]
                  else:
                      ret[idx_far] = val[pix_idx[idx_far], samp_idx[idx_far], -1]
  
              # d(nearest distance slice) < d < d(farthest distance slice)
              idx_btw = ~idx_near & ~idx_far & in_bounds_idx
              if np.any(idx_btw):
                  DM_ceil = self._DM_bin_edges[bin_idx_ceil[idx_btw]]
                  DM_floor = self._DM_bin_edges[bin_idx_ceil[idx_btw]-1]
                  a = (DM_ceil - dm[idx_btw]) / (DM_ceil - DM_floor)
                  if isinstance(samp_idx, slice):
                      ret[idx_btw] = (
                          (1.-a[:,None])
                          * val[pix_idx[idx_btw], samp_idx, bin_idx_ceil[idx_btw]]
                          + a[:,None]
                          * val[pix_idx[idx_btw], samp_idx, bin_idx_ceil[idx_btw]-1]
                      )
                  else:
                      ret[idx_btw] = (
                          (1.-a) * val[pix_idx[idx_btw], samp_idx[idx_btw], bin_idx_ceil[idx_btw]]
                          +    a * val[pix_idx[idx_btw], samp_idx[idx_btw], bin_idx_ceil[idx_btw]-1]
                      )
  
              # Flag: distance in reliable range?
              if return_flags:
                  dm_min = self._pixel_info['DM_reliable_min'][pix_idx]
                  dm_max = self._pixel_info['DM_reliable_max'][pix_idx]
                  flags['reliable_dist'] = (
                      (dm >= dm_min) &
                      (dm <= dm_max) &
                      np.isfinite(dm_min) &
                      np.isfinite(dm_max))
                  flags['reliable_dist'][~in_bounds_idx] = False
          else:   # No distances provided
              ret = val[pix_idx, samp_idx, :]   # Return all distances
              ret[~in_bounds_idx] = np.nan
  
              # Flag: reliable distance bounds
              if return_flags:
                  dm_min = self._pixel_info['DM_reliable_min'][pix_idx]
                  dm_max = self._pixel_info['DM_reliable_max'][pix_idx]
  
                  flags['min_reliable_distmod'] = dm_min
                  flags['max_reliable_distmod'] = dm_max
                  flags['min_reliable_distmod'][~in_bounds_idx] = np.nan
                  flags['max_reliable_distmod'][~in_bounds_idx] = np.nan
  
          # t4 = time.time()
  
          # Flag: convergence
          if return_flags:
              flags['converged'] = (
                  self._pixel_info['converged'][pix_idx].astype(np.bool))
              flags['converged'][~in_bounds_idx] = False
  
          # t5 = time.time()
  
          # Reduce the samples in the requested manner
          if mode == 'median':
              ret = np.median(ret, axis=1)
          elif mode == 'mean':
              ret = np.mean(ret, axis=1)
          elif mode == 'percentile':
              ret = np.nanpercentile(ret, pct, axis=1)
              if not scalar_pct:
                  # (percentile, pixel) -> (pixel, percentile)
                  # (pctile, pixel, distance) -> (pixel, distance, pctile)
                  ret = np.moveaxis(ret, 0, -1)
          elif mode == 'best':
              # Remove "samples" axis
              s = ret.shape
              ret.shape = s[:1] + s[2:]
          elif mode == 'samples':
              # Swap sample and distance axes to be consistent with other 3D dust
              # maps. The output shape will be (pixel, distance, sample).
              if not has_dist:
                  np.swapaxes(ret, 1, 2)
  
          # t6 = time.time()
          #
          # print('')
          # print('time inside bayestar.query: {:.4f} s'.format(t6-t0))
          # print('{: >7.4f} s : {: >6.4f} s : _find_data_idx'.format(t1-t0, t1-t0))
          # print('{: >7.4f} s : {: >6.4f} s : sample slice spec'.format(t2-t0, t2-t1))
          # print('{: >7.4f} s : {: >6.4f} s : create empty return flag array'.format(t3-t0, t3-t2))
          # print('{: >7.4f} s : {: >6.4f} s : extract results'.format(t4-t0, t4-t3))
          # print('{: >7.4f} s : {: >6.4f} s : convergence flag'.format(t5-t0, t5-t4))
          # print('{: >7.4f} s : {: >6.4f} s : reduce'.format(t6-t0, t6-t5))
          # print('')
  
          if return_flags:
              return ret, flags
  
          return ret