docs: def search(self, coords=None, distance_modulus=None, tolerance=1.e-2):
          """
          Organize a grid search over ROI target pixels and distance moduli in distance_modulus_array
          coords: (lon,lat)
          distance_modulus: scalar
          """
          nmoduli = len(self.distance_modulus_array)
          npixels    = len(self.roi.pixels_target)
          self.log_likelihood_sparse_array       = np.zeros([nmoduli, npixels])
          self.richness_sparse_array             = np.zeros([nmoduli, npixels])
          self.richness_lower_sparse_array       = np.zeros([nmoduli, npixels])
          self.richness_upper_sparse_array       = np.zeros([nmoduli, npixels])
          self.richness_upper_limit_sparse_array = np.zeros([nmoduli, npixels])
          self.stellar_mass_sparse_array         = np.zeros([nmoduli, npixels])
          self.fraction_observable_sparse_array  = np.zeros([nmoduli, npixels])
  
          # Specific pixel/distance_modulus
          coord_idx, distance_modulus_idx = None, None
          if coords is not None:
              # Match to nearest grid coordinate index
              coord_idx = self.roi.indexTarget(coords[0],coords[1])
          if distance_modulus is not None:
              # Match to nearest distance modulus index
              distance_modulus_idx=np.fabs(self.distance_modulus_array-distance_modulus).argmin()
  
          lon, lat = self.roi.pixels_target.lon, self.roi.pixels_target.lat
              
          logger.info('Looping over distance moduli in grid search ...')
          for ii, distance_modulus in enumerate(self.distance_modulus_array):
  
              # Specific pixel
              if distance_modulus_idx is not None:
                  if ii != distance_modulus_idx: continue
  
              logger.info('  (%-2i/%i) Distance Modulus=%.1f ...'%(ii+1,nmoduli,distance_modulus))
  
              # Set distance_modulus once to save time
              self.loglike.set_params(distance_modulus=distance_modulus)
  
              for jj in range(0, npixels):
                  # Specific pixel
                  if coord_idx is not None:
                      if jj != coord_idx: continue
  
                  # Set kernel location
                  self.loglike.set_params(lon=lon[jj],lat=lat[jj])
                  # Doesn't re-sync distance_modulus each time
                  self.loglike.sync_params()
                                           
                  args = (jj+1, npixels, self.loglike.source.lon, self.loglike.source.lat)
                  msg = '    (%-3i/%i) Candidate at (%.2f, %.2f) ... '%(args)
  
                  self.log_likelihood_sparse_array[ii][jj], self.richness_sparse_array[ii][jj], parabola = self.loglike.fit_richness()
                  self.stellar_mass_sparse_array[ii][jj] = self.stellar_mass_conversion * self.richness_sparse_array[ii][jj]
                  self.fraction_observable_sparse_array[ii][jj] = self.loglike.f
                  if self.config['scan']['full_pdf']:
                      #n_pdf_points = 100
                      #richness_range = parabola.profileUpperLimit(delta=25.) - self.richness_sparse_array[ii][jj]
                      #richness = np.linspace(max(0., self.richness_sparse_array[ii][jj] - richness_range),
                      #                          self.richness_sparse_array[ii][jj] + richness_range,
                      #                          n_pdf_points)
                      #if richness[0] > 0.:
                      #    richness = np.insert(richness, 0, 0.)
                      #    n_pdf_points += 1
                      # 
                      #log_likelihood = np.zeros(n_pdf_points)
                      #for kk in range(0, n_pdf_points):
                      #    log_likelihood[kk] = self.loglike.value(richness=richness[kk])
                      #parabola = ugali.utils.parabola.Parabola(richness, 2.*log_likelihood)
                      #self.richness_lower_sparse_array[ii][jj], self.richness_upper_sparse_array[ii][jj] = parabola.confidenceInterval(0.6827)
                      self.richness_lower_sparse_array[ii][jj], self.richness_upper_sparse_array[ii][jj] = self.loglike.richness_interval(0.6827)
                      
                      self.richness_upper_limit_sparse_array[ii][jj] = parabola.bayesianUpperLimit(0.95)
  
                      args = (
                          2. * self.log_likelihood_sparse_array[ii][jj],
                          self.stellar_mass_conversion*self.richness_sparse_array[ii][jj],
                          self.stellar_mass_conversion*self.richness_lower_sparse_array[ii][jj],
                          self.stellar_mass_conversion*self.richness_upper_sparse_array[ii][jj],
                          self.stellar_mass_conversion*self.richness_upper_limit_sparse_array[ii][jj]
                      )
                      msg += 'TS=%.1f, Stellar Mass=%.1f (%.1f -- %.1f @ 0.68 CL, < %.1f @ 0.95 CL)'%(args)
                  else:
                      args = (
                          2. * self.log_likelihood_sparse_array[ii][jj], 
                          self.stellar_mass_conversion * self.richness_sparse_array[ii][jj],
                          self.fraction_observable_sparse_array[ii][jj]
                      )
                      msg += 'TS=%.1f, Stellar Mass=%.1f, Fraction=%.2g'%(args)
                  logger.debug(msg)
                  
                  #if coords is not None and distance_modulus is not None:
                  #    results = [self.richness_sparse_array[ii][jj],
                  #               self.log_likelihood_sparse_array[ii][jj],
                  #               self.richness_lower_sparse_array[ii][jj],
                  #               self.richness_upper_sparse_array[ii][jj],
                  #               self.richness_upper_limit_sparse_array[ii][jj],
                  #               richness, log_likelihood, self.loglike.p, self.loglike.f]
                  #    return results
  
              jj_max = self.log_likelihood_sparse_array[ii].argmax()
              args = (
                  jj_max+1, npixels, lon[jj_max], lat[jj_max],
                  2. * self.log_likelihood_sparse_array[ii][jj_max], 
                  self.stellar_mass_conversion * self.richness_sparse_array[ii][jj_max]
              )
              msg = '  (%-3i/%i) Maximum at (%.2f, %.2f) ... TS=%.1f, Stellar Mass=%.1f'%(args)
              logger.info(msg)