Organize a grid search over ROI target pixels and distance moduli in distance_modulus_array
        coords: (lon,lat)
        distance_modulus: scalar

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)