def pixelizeCatalog(infiles, config, force=False):
    """
    Break catalog up into a set of healpix files.
    """
    nside_catalog = config['coords']['nside_catalog']
    nside_pixel = config['coords']['nside_pixel']
    outdir = mkdir(config['catalog']['dirname'])
    filenames = config.getFilenames()
    
    for ii,infile in enumerate(infiles):
        logger.info('(%i/%i) %s'%(ii+1, len(infiles), infile))
        f = pyfits.open(infile)
        data = f[1].data
        header = f[1].header
        logger.info("%i objects found"%len(data))
        if not len(data): continue
        glon,glat = cel2gal(data['RA'],data['DEC'])
        catalog_pix = ang2pix(nside_catalog,glon,glat,coord='GAL')
        pixel_pix = ang2pix(nside_pixel,glon,glat,coord='GAL')
        names = [n.upper() for n in data.columns.names]
        ra_idx = names.index('RA'); dec_idx = names.index('DEC')
        idx = ra_idx if ra_idx > dec_idx else dec_idx
        catalog_pix_name = 'PIX%i'%nside_catalog
        pixel_pix_name = 'PIX%i'%nside_pixel

        coldefs = pyfits.ColDefs(
            [pyfits.Column(name='GLON',format='1D',array=glon),
             pyfits.Column(name='GLAT',format='1D',array=glat),
             pyfits.Column(name=catalog_pix_name,format='1J',array=catalog_pix),
             pyfits.Column(name=pixel_pix_name  ,format='1J',array=pixel_pix)]
        )
        hdu = pyfits.new_table(data.columns[:idx+1]+coldefs+data.columns[idx+1:])
        table = hdu.data

        for pix in numpy.unique(catalog_pix):
            logger.debug("Processing pixel %s"%pix)
            outfile = filenames.data['catalog'][pix]
            if not os.path.exists(outfile):
                logger.debug("Creating %s"%outfile)
                names = [n.upper() for n in table.columns.names]
                formats = table.columns.formats
                columns = [pyfits.Column(n,f) for n,f in zip(names,formats)]
                out = pyfits.HDUList([pyfits.PrimaryHDU(),pyfits.new_table(columns)])
                out[1].header['NSIDE'] = nside_catalog
                out[1].header['PIX'] = pix
                out.writeto(outfile)
            hdulist = pyfits.open(outfile,mode='update')
            t1 = hdulist[1].data
            # Could we speed up with sorting and indexing?
            t2 = table[ table[catalog_pix_name] == pix ]
            nrows1 = t1.shape[0]
            nrows2 = t2.shape[0]
            nrows = nrows1 + nrows2
            out = pyfits.new_table(t1.columns, nrows=nrows)
            for name in t1.columns.names:
                out.data.field(name)[nrows1:]=t2.field(name)
            hdulist[1] = out
            logger.debug("Writing %s"%outfile)
            hdulist.flush()
            hdulist.close()

def run(self):
    if 'pixelize' in self.opts.run:
        # Pixelize the raw catalog data
        logger.info("Running 'pixelize'...")
        rawdir = self.config['data']['dirname']
        rawfiles = sorted(glob.glob(os.path.join(rawdir,'*.fits')))
        x = ugali.preprocess.pixelize.pixelizeCatalog(rawfiles,self.config)
    if 'density' in self.opts.run:
        # Calculate magnitude limits
        logger.info("Running 'density'...")
        x = ugali.preprocess.pixelize.pixelizeDensity(self.config,nside=2**9,force=self.opts.force)
    if 'maglims' in self.opts.run:
        # Calculate magnitude limits
        logger.info("Running 'maglims'...")
        maglims = ugali.preprocess.maglims.Maglims(self.config)
        x = maglims.run(force=self.opts.force)
    if 'simple' in self.opts.run:
        # Calculate simple magnitude limits
        logger.info("Running 'simple'...")
        #ugali.preprocess.maglims.simple_maglims(self.config,force=self.opts.force)
        maglims = ugali.preprocess.maglims.Maglims(self.config)
        x = maglims.run(simple=True,force=self.opts.force)
    if 'split' in self.opts.run:
        logger.info("Running 'split'...")
        ugali.preprocess.maglims.simple_split(self.config,'split',force=self.opts.force)

    def download(self, pixel, outdir=None, force=False):
        import pyfits 

        if outdir is None: outdir = './'
        else:              mkdir(outdir)
        sqldir = mkdir(os.path.join(outdir,'sql'))
        self._setup_desdbi()

        basename = self.basename + "_%04d"%pixel['name']
        sqlname = os.path.join(sqldir,basename+'.sql')
        taskname = basename
        outfile = os.path.join(outdir,basename+".fits")
        # ADW: There should be a 'force' option here
        if os.path.exists(outfile) and not force:
            logger.warning("Found %s; skipping..."%(outfile))
            return

        logger.info("\nDownloading pixel: %(name)i (ra=%(ra_min)g:%(ra_max)g,dec=%(dec_min)g:%(dec_max)g)"%(pixel))
        logger.info("Working on "+sqlname)
         
        self.generate_query(pixel['ra_min'],pixel['ra_max'],pixel['dec_min'],pixel['dec_max'],sqlname,outfile)
        ret = self.query(self.release,taskname,sqlname)
        if ret != 0:
            msg = "Download failed to complete."
            raise Exception(msg)
        return outfile

    def observableFractionMMD(self, mask, distance_modulus, mass_min=0.1):
        # This can be done faster...
        logger.info('Calculating observable fraction from MMD')

        mmd = self.signalMMD(mask,distance_modulus)
        obs_frac = mmd.sum(axis=-1).sum(axis=-1)[mask.mask_roi_digi[mask.roi.pixel_interior_cut]]
        return obs_frac

def pixelizeDensity(config, nside=None, force=False):
    if nside is None: 
        nside = config['coords']['nside_likelihood']
    coordsys = config['coords']['coordsys'].upper()
    filenames = config.getFilenames()
    infiles = filenames[~filenames['catalog'].mask]
    lon_field = config['catalog']['lon_field'].upper()
    lat_field = config['catalog']['lat_field'].upper()

    for ii,f in enumerate(infiles.data):
        infile = f['catalog']
        pix = f['pix']
        logger.info('(%i/%i) %s'%(ii+1, len(infiles), infile))

        outfile = config['data']['density']%pix
        if os.path.exists(outfile) and not force: 
            logger.info("Found %s; skipping..."%outfile)
            continue
            
        outdir = mkdir(os.path.dirname(outfile))
        pixels, density = stellarDensity(infile,nside,
                                         lon_field=lon_field,lat_field=lat_field)

        data = dict(PIXEL=pixels,DENSITY=density)
        healpix.write_partial_map(outfile,data,nside=nside,coord=coordsys[0])

    def pdf_mmd(self, lon, lat, mag_1, mag_2, distance_modulus, mask, delta_mag=0.03, steps=1000):
        """
        Ok, now here comes the beauty of having the signal MMD.
        """
        logger.info('Running MMD pdf')
 
        roi = mask.roi
        mmd = self.signalMMD(mask,distance_modulus,delta_mag=delta_mag,mass_steps=steps)
        
        # This is fragile, store this information somewhere else...
        nedges = np.rint((roi.bins_mag[-1]-roi.bins_mag[0])/delta_mag)+1
        edges_mag,delta_mag = np.linspace(roi.bins_mag[0],roi.bins_mag[-1],nedges,retstep=True)
                                    
        idx_mag_1 = np.searchsorted(edges_mag,mag_1)
        idx_mag_2 = np.searchsorted(edges_mag,mag_2)
 
        if np.any(idx_mag_1 > nedges) or np.any(idx_mag_1 == 0):
            msg = "Magnitude out of range..."
            raise Exception(msg)
        if np.any(idx_mag_2 > nedges) or np.any(idx_mag_2 == 0):
            msg = "Magnitude out of range..."
            raise Exception(msg)
 
        idx = mask.roi.indexROI(lon,lat)
        u_color = mmd[(mask.mask_roi_digi[idx],idx_mag_1,idx_mag_2)]
 
        # Remove the bin size to convert the pdf to units of mag^-2
        u_color /= delta_mag**2
 
        return u_color

    def _setup_cmd(self,mode='cloud-in-cells'):
        """
        The purpose here is to create a more finely binned
        background CMD to sample from.
        """
        # Only setup once...
        if hasattr(self,'bkg_lambda'): return

        logger.info("Setup color...")
        # In the limit theta->0: 2*pi*(1-cos(theta)) -> pi*theta**2
        # (Remember to convert from sr to deg^2) 
        #solid_angle_roi = sr2deg(2*np.pi*(1-np.cos(np.radians(self.roi_radius))))
        solid_angle_roi = self.roi.area_pixel*len(self.roi.pixels)

        # Large CMD bins cause problems when simulating
        config = Config(self.config) 
        config['color']['n_bins'] *= 5 #10
        config['mag']['n_bins']   *= 1 #2
        #config['mask']['minimum_solid_angle'] = 0
        roi = ugali.analysis.loglike.createROI(config,self.roi.lon,self.roi.lat)
        mask = ugali.analysis.loglike.createMask(config,roi)

        self.bkg_centers_color  = roi.centers_color
        self.bkg_centers_mag    = roi.centers_mag

        # Background CMD has units: [objs / deg^2 / mag^2]
        cmd_background = mask.backgroundCMD(self.catalog,mode)
        
        self.bkg_lambda=cmd_background*solid_angle_roi*roi.delta_color*roi.delta_mag
        np.sum(self.bkg_lambda)

        # Clean up 
        del config, roi, mask

 def query(self,dbase,task,query):
     logger.info("Running query...")
     cmd = "java -jar casjobs.jar run -t %s -n %s -f %s" % (dbase,task,query)
     logger.info(cmd)
     ret = subprocess.check_output(cmd,shell=True,stderr=subprocess.STDOUT) 
     if 'ERROR:' in ret:
         raise subprocess.CalledProcessError(1,cmd,ret)
     return ret

 def writeLabels(self,filename=None):
     if filename is None: filename = self.labelfile
     # Converting to float is a waste of memory...
     # This should be much more robustly done in writeSparseHealpixMap
     data_dict = {'LABEL':self.labels.astype(float)}
     logger.info("Writing %s..."%filename)
     ugali.utils.skymap.writeSparseHealpixMap(self.pixels,data_dict,self.nside,filename,
                                              distance_modulus_array=self.distances)

 def print_info(self,age,metallicity):
     params = dict(age=age,z=metallicity)
     params['name'] = self.__class__.__name__
     params['survey'] = self.survey
     params['feh'] = self.isochrone.z2feh(metallicity)
     msg = 'Downloading: %(name)s (survey=%(survey)s, age=%(age).1fGyr, Z=%(z).5f, Fe/H=%(feh).3f)'%params
     logger.info(msg)
     return msg

    def write(self, outfile):
        """
        Save the likelihood results as a sparse HEALPix map.
        """
        data = odict()
        data['PIXEL']=self.roi.pixels_target
        # Full data output (too large for survey)
        if self.config['scan']['full_pdf']:
            data['LOG_LIKELIHOOD']=self.log_likelihood_sparse_array.T
            data['RICHNESS']=self.richness_sparse_array.T
            data['RICHNESS_LOWER']=self.richness_lower_sparse_array.T
            data['RICHNESS_UPPER']=self.richness_upper_sparse_array.T
            data['RICHNESS_LIMIT']=self.richness_upper_limit_sparse_array.T
            #data['STELLAR_MASS']=self.stellar_mass_sparse_array.T
            data['FRACTION_OBSERVABLE']=self.fraction_observable_sparse_array.T
        else:
            data['LOG_LIKELIHOOD']=self.log_likelihood_sparse_array.T
            data['RICHNESS']=self.richness_sparse_array.T
            data['FRACTION_OBSERVABLE']=self.fraction_observable_sparse_array.T

        # Convert to 32bit float
        for k in list(data.keys())[1:]:
            data[k] = data[k].astype('f4',copy=False)
            
        # Stellar mass can be calculated from STELLAR * RICHNESS
        header = odict()
        header['STELLAR']=round(self.stellar_mass_conversion,8)
        header['LKDNSIDE']=self.config['coords']['nside_likelihood']
        header['LKDPIX']=ang2pix(self.config['coords']['nside_likelihood'],
                                 self.roi.lon,self.roi.lat)
        header['NROI']=self.roi.inROI(self.loglike.catalog_roi.lon,
                                      self.loglike.catalog_roi.lat).sum()
        header['NANNULUS']=self.roi.inAnnulus(self.loglike.catalog_roi.lon,
                                              self.loglike.catalog_roi.lat).sum()
        header['NINSIDE']=self.roi.inInterior(self.loglike.catalog_roi.lon,
                                              self.loglike.catalog_roi.lat).sum()
        header['NTARGET']=self.roi.inTarget(self.loglike.catalog_roi.lon,
                                            self.loglike.catalog_roi.lat).sum()

        # Flatten if there is only a single distance modulus
        # ADW: Is this really what we want to do?
        if len(self.distance_modulus_array) == 1:
            for key in data:
                data[key] = data[key].flatten()

        logger.info("Writing %s..."%outfile)
        write_partial_map(outfile,data,
                          nside=self.config['coords']['nside_pixel'],
                          header=header,
                          clobber=True
                          )
        
        fitsio.write(outfile,
                     dict(DISTANCE_MODULUS=self.distance_modulus_array.astype('f4',copy=False)),
                     extname='DISTANCE_MODULUS',
                     clobber=False)

 def throttle(self,max_jobs=None,sleep=60):
     if max_jobs is None: max_jobs = self.max_jobs
     if max_jobs is None: return
     while True:
         njobs = self.njobs()
         if njobs < max_jobs:
             return
         else:
             logger.info('%i jobs already in queue, waiting...'%(njobs))
             time.sleep(sleep)

    def writeCandidates(self,filename=None):
        if filename is None: filename = self.candfile

        threshold = self.config['search']['cand_threshold']
        select  = (self.assocs['CUT']==0)
        select &= (self.assocs['TS']>threshold)
        #select &= (self.assocs['ASSOC2']=='')

        self.candidates = self.assocs[select]
        logger.info("Writing %s..."%filename)
        fitsio.write(filename,self.candidates,clobber=True)

def do_membership(args):
    """ Write the membership output file """
    config,name,label,coord = args

    filenames = make_filenames(config,label)
    srcfile = filenames['srcfile']
    memfile = filenames['memfile']

    logger.info("Writing %s..."%memfile)
    from ugali.analysis.loglike import write_membership
    write_membership(memfile,config,srcfile,section='source')

    def labelHealpix(pixels, values, nside, threshold=0, xsize=1000):
        """
        Label contiguous regions of a (sparse) HEALPix map. Works by mapping 
        HEALPix array to a Mollweide projection and applying scipy.ndimage.label
     
        Assumes non-nested HEALPix map.
        
        Parameters:
        pixels    : Pixel values associated to (sparse) HEALPix array
        values    : (Sparse) HEALPix array of data values
        nside     : HEALPix dimensionality
        threshold : Threshold value for object detection
        xsize     : Size of Mollweide projection
        
        Returns:
        labels, nlabels
        """
        proj = healpy.projector.MollweideProj(xsize=xsize)
        vec = healpy.pix2vec(nside,pixels)
        xy = proj.vec2xy(vec)
        ij = proj.xy2ij(xy)
        xx,yy = proj.ij2xy()
     
        # Convert to Mollweide
        searchims = []
        if values.ndim < 2: iterate = [values]
        else:               iterate = values.T
        for i,value in enumerate(iterate):
            logger.debug("Labeling slice %i...")
            searchim = numpy.zeros(xx.shape,dtype=bool)
            select = (value > threshold)
            yidx = ij[0][select]; xidx = ij[1][select]
            searchim[yidx,xidx] |= True
            searchims.append( searchim )
        searchims = numpy.array(searchims)

        # Full binary structure
        s = ndimage.generate_binary_structure(searchims.ndim,searchims.ndim)
     
        ### # Dilate in the z-direction
        logger.info("  Dilating image...")
        searchims = ndimage.binary_dilation(searchims,s,1)
        
        # Do the labeling
        logger.info("  Labeling image...")
        labels,nlabels = ndimage.label(searchims,structure=s)

        # Convert back to healpix
        pix_labels = labels[:,ij[0],ij[1]].T
        pix_labels = pix_labels.reshape(values.shape)
        pix_labels *= (values > threshold) # re-trim

        return pix_labels, nlabels

 def sky(self,lon=None,lat=None,size=1):
     logger.info("Generating %i random points..."%size)
     # Random longitue and latitude
     lon,lat = ugali.utils.stats.sky(lon,lat,size=10*size)
     # Random healpix coordinates inside footprint
     nside_pixel = self.config['coords']['nside_pixel']
     pixels = ang2pix(nside_pixel,lon,lat)
     if np.unique(pixels).size > 1:
         inside = ugali.utils.skymap.inFootprint(self.config,pixels,nside=nside_pixel)
     else:
         inside = np.ones(len(pixels),dtype=bool)
     return lon[inside][:size],lat[inside][:size]

    def writeLabels(self,filename=None):
        if filename is None: filename = self.labelfile

        # ADW: Is it necessary to convert labels?
        data_dict = {'PIXEL':self.pixels,
                     'LABEL':self.labels.astype(float,copy=False)}

        logger.info("Writing %s..."%filename)
        healpix.write_partial_map(filename,data_dict,self.nside)
        fitsio.write(filename,
                     {'DISTANCE_MODULUS':self.distances.astype('f4',copy=False)},
                     extname='DISTANCE_MODULUS',
                     clobber=False)

    def _applySelection(self,selection=None):
        # ADW: This is a hack (eval is unsafe!)
        if selection is None:
            selection = self.config['catalog'].get('selection')

        if not selection: 
            return
        elif 'self.data' not in selection:
            msg = "Selection does not contain 'data'"
            raise Exception(msg)
        else:
            logger.info('Evaluating selection: \n"%s"'%selection)
            sel = eval(selection)
            self.data = self.data[sel]

    def spatialBin(self, roi):
        """
        Calculate indices of ROI pixels corresponding to object locations.
        """
        if hasattr(self,'pixel_roi_index') and hasattr(self,'pixel'): 
            logger.warning('Catalog alread spatially binned')
            return

        # ADW: Not safe to set index = -1 (since it will access last entry); 
        # np.inf would be better...
        self.pixel = ang2pix(self.config['coords']['nside_pixel'],self.lon,self.lat)
        self.pixel_roi_index = roi.indexROI(self.lon,self.lat)

        logger.info("Found %i objects outside ROI"%(self.pixel_roi_index < 0).sum())

    def satellite(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwargs):
        """
        Create a simulated satellite. Returns a catalog object.
        """
        if seed is not None: np.random.seed(seed)

        isochrone = kwargs.pop('isochrone',self.isochrone)
        kernel    = kwargs.pop('kernel',self.kernel)

        for k,v in kwargs.items():
            if k in kernel.params.keys(): setattr(kernel,k,v)

        mag_1, mag_2 = isochrone.simulate(stellar_mass, distance_modulus)
        lon, lat     = kernel.simulate(len(mag_1))
 
        logger.info("Simulating %i satellite stars..."%len(mag_1))
        pix = ang2pix(self.config['coords']['nside_pixel'], lon, lat)

        # There is probably a better way to do this step without creating the full HEALPix map
        mask = -1. * numpy.ones(healpy.nside2npix(self.config['coords']['nside_pixel']))
        mask[self.roi.pixels] = self.mask.mask_1.mask_roi_sparse
        mag_lim_1 = mask[pix]
        mask = -1. * numpy.ones(healpy.nside2npix(self.config['coords']['nside_pixel']))
        mask[self.roi.pixels] = self.mask.mask_2.mask_roi_sparse
        mag_lim_2 = mask[pix]

        mag_err_1 = self.photo_err_1(mag_lim_1 - mag_1)
        mag_err_2 = self.photo_err_2(mag_lim_2 - mag_2)

        # Randomize magnitudes by their errors
        mag_obs_1 = mag_1+numpy.random.normal(size=len(mag_1))*mag_err_1
        mag_obs_2 = mag_2+numpy.random.normal(size=len(mag_2))*mag_err_2
        #mag_obs_1 = mag_1
        #mag_obs_2 = mag_2

        #select = numpy.logical_and(mag_obs_1 < mag_lim_1, mag_obs_2 < mag_lim_2)
        select = (mag_lim_1>mag_obs_1)&(mag_lim_2>mag_obs_2)

        # Make sure objects lie within the original cmd (should also be done later...)
        #select &= (ugali.utils.binning.take2D(self.mask.solid_angle_cmd, mag_obs_1 - mag_obs_2, mag_obs_1,self.roi.bins_color, self.roi.bins_mag) > 0)

        #return mag_1_obs[cut], mag_2_obs[cut], lon[cut], lat[cut]
        logger.info("Clipping %i simulated satellite stars..."%(~select).sum())
        mc_source_id = mc_source_id * numpy.ones(len(mag_1))
        
        hdu = ugali.observation.catalog.makeHDU(self.config,mag_obs_1[select],mag_err_1[select],
                                                mag_obs_2[select],mag_err_2[select], 
                                                lon[select],lat[select],mc_source_id[select])
        catalog = ugali.observation.catalog.Catalog(self.config, data=hdu.data)
        return catalog