diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index fa862a1..c8b745b 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -7,9 +7,9 @@ name: build
 
 on:
   push:
-    branches: [ master, actions ]
   pull_request:
-    branches: [ master ]
+  schedule:
+  - cron: '0 0 1 * *'
 
 jobs:
   build:
@@ -38,19 +38,17 @@ jobs:
         source activate env
         export UGALIDIR="$HOME/.ugali"
         python setup.py -q install --isochrones --catalogs
-    - name: Install test data
-      run: |
-        wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-test-data.tar.gz -O ugali-test-data.tar.gz
-        tar -xzf ugali-test-data.tar.gz
     - name: Lint with flake8
-      if: ${{ false }}
+      if: ${{ true }}
       run: |
         source activate env
         conda install -q flake8 -c conda-forge
         # stop the build if there are Python syntax errors or undefined names
         flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
-        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    - name: Install test data
+      run: |
+        wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-test-data.tar.gz -O ugali-test-data.tar.gz
+        tar -xzf ugali-test-data.tar.gz
     - name: Test with nose
       run: |
         source activate env
diff --git a/ugali/analysis/prior.py b/ugali/analysis/prior.py
index 58ac621..204b93f 100644
--- a/ugali/analysis/prior.py
+++ b/ugali/analysis/prior.py
@@ -1,4 +1,9 @@
 #!/usr/bin/env python
+"""
+Incomplete module for dealing with priors...
+"""
+
+import scipy.stats
 
 class Prior(object):
     def __call__(self, value):
diff --git a/ugali/analysis/results.py b/ugali/analysis/results.py
index c98796c..99877d6 100644
--- a/ugali/analysis/results.py
+++ b/ugali/analysis/results.py
@@ -18,7 +18,8 @@
 
 import ugali.utils.stats
 from ugali.utils.stats import Samples
-from ugali.utils.projector import dist2mod,mod2dist,gal2cel,gal2cel_angle
+from ugali.utils.projector import dist2mod,mod2dist
+from ugali.utils.projector import cel2gal,gal2cel,gal2cel_angle
 from ugali.utils.projector import ang2const, ang2iau
 from ugali.utils.config import Config
 from ugali.utils.logger import logger
diff --git a/ugali/analysis/scan.py b/ugali/analysis/scan.py
index 93c1ea3..e3dca85 100755
--- a/ugali/analysis/scan.py
+++ b/ugali/analysis/scan.py
@@ -130,10 +130,12 @@ def search(self, coords=None, distance_modulus=None, extension=None, tolerance=1
         self.stellar_mass_array         = np.zeros([nmoduli,npixels],dtype='f4')
         self.fraction_observable_array  = np.zeros([nmoduli,npixels],dtype='f4')
         self.extension_fit_array        = np.zeros([nmoduli,npixels],dtype='f4')
-        # DEPRECATED: ADW 2019-04-27
-        self.richness_lower_array       = np.zeros([nmoduli,npixels],dtype='f4')
-        self.richness_upper_array       = np.zeros([nmoduli,npixels],dtype='f4')
-        self.richness_ulimit_array      = np.zeros([nmoduli,npixels],dtype='f4')
+        if self.config['scan']['full_pdf']:
+            # DEPRECATED: ADW 2019-04-27
+            DeprecationWarning("'full_pdf' is deprecated.")
+            self.richness_lower_array       = np.zeros([nmoduli,npixels],dtype='f4')
+            self.richness_upper_array       = np.zeros([nmoduli,npixels],dtype='f4')
+            self.richness_ulimit_array      = np.zeros([nmoduli,npixels],dtype='f4')
 
         # Specific pixel/distance_modulus
         coord_idx, distance_modulus_idx, extension_idx = None, None, None
@@ -315,6 +317,7 @@ def write(self, outfile):
         data['PIXEL']=self.roi.pixels_target
         # Full data output (too large for survey)
         if self.config['scan']['full_pdf']:
+            DeprecationWarning("'full_pdf' is deprecated.")
             data['LOG_LIKELIHOOD']=self.loglike_array.T
             data['RICHNESS']=self.richness_array.T
             data['RICHNESS_LOWER']=self.richness_lower_array.T
diff --git a/ugali/candidate/associate.py b/ugali/candidate/associate.py
index 4ca91ce..1667ee8 100644
--- a/ugali/candidate/associate.py
+++ b/ugali/candidate/associate.py
@@ -54,7 +54,7 @@
 #            glon, glat = ugali.utils.projector.celToGal(lon,lat)
 #        else:
 #            glon,glat = lon, lat
-#        return ugali.utils.projector.match(glon,glat,self.data['glon'],self.data['glat'],tol)
+#        return ugali.utils.projector.match(glon,glat,self['glon'],self['glat'],tol)
 
 
 
@@ -452,7 +452,7 @@ def catalogFactory(name, **kwargs):
     catalogs = odict(inspect.getmembers(sys.modules[__name__], fn))
 
     if name not in list(catalogs.keys()):
-        msg = "%s not found in catalogs:\n %s"%(name,list(kernels.keys()))
+        msg = "%s not found in catalogs:\n %s"%(name,list(catalogs.keys()))
         logger.error(msg)
         msg = "Unrecognized catalog: %s"%name
         raise Exception(msg)
diff --git a/ugali/isochrone/dartmouth.py b/ugali/isochrone/dartmouth.py
index 3e0a2d1..7cf713b 100644
--- a/ugali/isochrone/dartmouth.py
+++ b/ugali/isochrone/dartmouth.py
@@ -165,7 +165,7 @@ def _parse(self,filename):
         try:
             columns = self.columns[self.survey.lower()]
         except KeyError as e:
-            logger.warning('Unrecognized survey: %s'%(survey))
+            logger.warning('Unrecognized survey: %s'%(self.survey))
             raise(e)
 
         kwargs = dict(comments='#',usecols=list(columns.keys()),dtype=list(columns.values()))
diff --git a/ugali/isochrone/empirical.py b/ugali/isochrone/empirical.py
index 06b316f..b309a13 100644
--- a/ugali/isochrone/empirical.py
+++ b/ugali/isochrone/empirical.py
@@ -3,7 +3,10 @@
 Generic python script.
 """
 import os
+from collections import OrderedDict as odict
 
+from ugali.analysis.model import Model, Parameter
+from ugali.utils.shell import mkdir, get_ugali_dir, get_iso_dir
 from ugali.isochrone.parsec import PadovaIsochrone
 
 class EmpiricalPadova(PadovaIsochrone):
diff --git a/ugali/isochrone/mesa.py b/ugali/isochrone/mesa.py
index 70df82b..28ca103 100644
--- a/ugali/isochrone/mesa.py
+++ b/ugali/isochrone/mesa.py
@@ -134,7 +134,7 @@ def _parse(self,filename):
         try:
             columns = self.columns[self.survey.lower()]
         except KeyError as e:
-            logger.warning('Unrecognized survey: %s'%(survey))
+            logger.warning('Unrecognized survey: %s'%(self.survey))
             raise(e)
 
         kwargs = dict(comments='#',usecols=list(columns.keys()),dtype=list(columns.values()))
diff --git a/ugali/isochrone/model.py b/ugali/isochrone/model.py
index 8f592dc..535eb3f 100644
--- a/ugali/isochrone/model.py
+++ b/ugali/isochrone/model.py
@@ -322,27 +322,6 @@ def stellar_luminosity(self, steps=10000):
 
         return np.sum(luminosity * d_log_mass * self.imf.pdf(mass, log_mode=True))
 
-    def stellar_luminosity2(self, steps=10000):
-        """
-        DEPRECATED: ADW 2017-09-20
-
-        Compute the stellar luminosity (L_Sol; average per star).
-        Uses "sample" to generate mass sample and pdf.  The range of
-        integration only covers the input isochrone data (no
-        extrapolation used), but this seems like a sub-percent effect
-        if the isochrone goes to 0.15 Msun for the old and metal-poor
-        stellar populations of interest.
-
-        Note that the stellar luminosity is very sensitive to the
-        post-AGB population.
-        """
-        msg = "'%s.stellar_luminosity2': ADW 2017-09-20"%self.__class__.__name__
-        DeprecationWarning(msg)
-        mass_init, mass_pdf, mass_act, mag_1, mag_2 = self.sample(mass_steps=steps)
-        luminosity_interpolation = scipy.interpolate.interp1d(self.mass_init, self.luminosity,fill_value=0,bounds_error=False)
-        luminosity = luminosity_interpolation(mass_init)
-        return np.sum(luminosity * mass_pdf)
-
     # ADW: For temporary backward compatibility
     stellarMass = stellar_mass
     stellarLuminosity = stellar_luminosity
@@ -354,7 +333,6 @@ def absolute_magnitude(self, richness=1, steps=1e4):
         g,r -> V transform equations from Jester 2005
         [astro-ph/0506022].
 
-
         TODO: ADW If richness not specified, should use self.richness
 
         Parameters:
diff --git a/ugali/isochrone/padova.py b/ugali/isochrone/padova.py
index 9c919ad..b6236c3 100644
--- a/ugali/isochrone/padova.py
+++ b/ugali/isochrone/padova.py
@@ -1,7 +1,20 @@
 #!/usr/bin/env python
 """
 Older (pre-PARSEC) Padova isochrones.
+
+Not verified to work...
 """
+import os.path
+from collections import OrderedDict as odict
+
+import numpy as np
+
+from ugali.analysis.model import Model, Parameter
+from ugali.utils.shell import mkdir, get_ugali_dir, get_iso_dir
+from ugali.isochrone.model import Isochrone
+from ugali.isochrone.parsec import Marigo2017 as Padova
+from ugali.isochrone.parsec import defaults_27
+from ugali.utils.logger import logger
 
 class Girardi2002(Padova):
     defaults = dict(defaults_27)
@@ -19,7 +32,7 @@ class Girardi2010b(Padova):
     defaults = dict(defaults_27)
     defaults['isoc_kind'] = 'gi10b'
 
-class Girardi2002(PadovaIsochrone):
+class Girardi2002(Padova):
     #_dirname = '/u/ki/kadrlica/des/isochrones/v5/'
     _dirname =  os.path.join(get_iso_dir(),'{survey}','girardi2002')
     # For use with Marigo et al. (2008) and earlier use Anders & Grevesse 1989
@@ -57,7 +70,7 @@ def _parse(self,filename):
         try:
             columns = self.columns[self.survey.lower()]
         except KeyError as e:
-            logger.warning('did not recognize survey %s'%(survey))
+            logger.warning('did not recognize survey %s'%(self.survey))
             raise(e)
 
         kwargs = dict(delimiter='\t',usecols=list(columns.keys()),dtype=list(columns.values()))
diff --git a/ugali/isochrone/parsec.py b/ugali/isochrone/parsec.py
index 5157fec..5f784e3 100644
--- a/ugali/isochrone/parsec.py
+++ b/ugali/isochrone/parsec.py
@@ -370,7 +370,7 @@ def _parse(self,filename):
         try:
             columns = self.columns[self.survey.lower()]
         except KeyError as e:
-            logger.warning('Unrecognized survey: %s'%(survey))
+            logger.warning('Unrecognized survey: %s'%(self.survey))
             raise(e)
 
         # delimiter='\t' is used to be compatible with OldPadova...
@@ -455,7 +455,7 @@ class Marigo2017(ParsecIsochrone):
                 (28,('Y',float)),
                 ]),
         )
-    
+
     def _parse(self,filename):
         """Reads an isochrone file in the Padova (Marigo et al. 2017)
         format. Creates arrays with the initial stellar mass and
@@ -464,7 +464,7 @@ def _parse(self,filename):
         try:
             columns = self.columns[self.survey.lower()]
         except KeyError as e:
-            logger.warning('Unrecognized survey: %s'%(survey))
+            logger.warning('Unrecognized survey: %s'%(self.survey))
             raise(e)
 
         kwargs = dict(usecols=list(columns.keys()),dtype=list(columns.values()))
diff --git a/ugali/observation/mask.py b/ugali/observation/mask.py
index 5ad6242..c59146e 100644
--- a/ugali/observation/mask.py
+++ b/ugali/observation/mask.py
@@ -13,7 +13,6 @@
 import healpy as hp
 import scipy.signal
 
-#import ugali.utils.plotting
 import ugali.utils.binning
 import ugali.utils.skymap
 
@@ -403,51 +402,19 @@ def photo_err_2(delta):
 
         return self.photo_err_1, self.photo_err_2
 
-    def plotSolidAngleCMD(self):
-        """
-        Solid angle within the mask as a function of color and magnitude.
-        """
-        msg = "'%s.plotSolidAngleCMD': ADW 2018-05-05"%self.__class__.__name__
-        DeprecationWarning(msg)
-
-        import ugali.utils.plotting        
-        ugali.utils.plotting.twoDimensionalHistogram('mask', 'color', 'magnitude',
-                                                     self.solid_angle_cmd,
-                                                     self.roi.bins_color,
-                                                     self.roi.bins_mag,
-                                                     lim_x = [self.roi.bins_color[0],
-                                                              self.roi.bins_color[-1]],
-                                                     lim_y = [self.roi.bins_mag[-1],
-                                                              self.roi.bins_mag[0]])
-
-    def plotSolidAngleMMD(self):
-        """
-        Solid angle within the mask as a function of color and magnitude.
-        """
-        msg = "'%s.plotSolidAngleMMD': ADW 2018-05-05"%self.__class__.__name__
-        DeprecationWarning(msg)
-
-        import ugali.utils.plotting        
-
-        ugali.utils.plotting.twoDimensionalHistogram('mask', 'mag_2', 'mag_1',
-                                                     self.solid_angle_mmd,
-                                                     self.roi.bins_mag,
-                                                     self.roi.bins_mag,
-                                                     lim_x = [self.roi.bins_mag[0],
-                                                              self.roi.bins_mag[-1]],
-                                                     lim_y = [self.roi.bins_mag[-1],
-                                                              self.roi.bins_mag[0]])
-
-
-
     def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None):
         """
         Generate an empirical background model in magnitude-magnitude space.
         
-        INPUTS:
-            catalog: Catalog object
-        OUTPUTS:
-            background
+        Parameters
+        ----------
+        catalog: catalog object
+        method:  method for estimated MMD
+        weights: weights assigned to each catalog object
+
+        Returns
+        -------
+        background: estimate of background magnitude-magnitude distribution
         """
 
         # Select objects in annulus
@@ -474,10 +441,8 @@ def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None):
                                                               [self.roi.bins_mag,self.roi.bins_mag],
                                                               weights=number_density)[0]
         elif method == 'bootstrap':
-            # Not implemented
+            # Not implemented; see catalog.bootstrap
             raise ValueError("Bootstrap method not implemented")
-            mag_1 + (mag_1_err * np.random.normal(0, 1., len(mag_1)))
-            mag_2 + (mag_2_err * np.random.normal(0, 1., len(mag_2)))
 
         elif method == 'histogram':
             # Apply raw histogram
@@ -547,10 +512,15 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None):
         """
         Generate an empirical background model in color-magnitude space.
         
-        INPUTS:
-            catalog: Catalog object
-        OUTPUTS:
-            background
+        Parameters
+        ----------
+        catalog: catalog object
+        method:  method for estimated MMD
+        weights: weights assigned to each catalog object
+
+        Returns
+        -------
+        background: estimate of background color-magnitude distribution
         """
 
         # Select objects in annulus
@@ -577,13 +547,8 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None):
                                                               [self.roi.bins_color,self.roi.bins_mag],
                                                               weights=number_density)[0]
         elif mode == 'bootstrap':
-            # Not implemented
+            # Not implemented; see catalog.bootstrap
             raise ValueError("Bootstrap mode not implemented")
-            mag_1_array = catalog.mag_1
-            mag_2_array = catalog.mag_2
-
-            catalog.mag_1 + (catalog.mag_1_err * np.random.normal(0, 1., len(catalog.mag_1)))
-            catalog.mag_2 + (catalog.mag_2_err * np.random.normal(0, 1., len(catalog.mag_2)))
 
         elif mode == 'histogram':
             # Apply raw histogram
@@ -847,24 +812,6 @@ def depth(self, lon, lat):
         """
         pass
 
-    def plot(self):
-        """
-        Plot the magnitude depth.
-        """
-        msg = "'%s.plot': ADW 2018-05-05"%self.__class__.__name__
-        DeprecationWarning(msg)
-
-        import ugali.utils.plotting
-
-        mask = hp.UNSEEN * np.ones(hp.nside2npix(self.nside))
-        mask[self.roi.pixels] = self.mask_roi_sparse
-        mask[mask == 0.] = hp.UNSEEN
-        ugali.utils.plotting.zoomedHealpixMap('Completeness Depth',
-                                              mask,
-                                              self.roi.lon, self.roi.lat,
-                                              self.roi.config.params['coords']['roi_radius'])
-
-
 class CoverageBand(object):
     """
     Map of coverage fraction for a single observing band.
@@ -923,106 +870,3 @@ def __init__(self, maglim, roi):
         self.mask_roi_sparse = mask[self.roi.pixels] 
 
 ############################################################
-def simpleMask(config):
-
-    #params = ugali.utils.(config, kwargs)
-
-    roi = ugali.observation.roi.ROI(config)
-
-    # De-project the bin centers to get magnitude depths
-
-    mesh_x, mesh_y = np.meshgrid(roi.centers_x, roi.centers_y)
-    r = np.sqrt(mesh_x**2 + mesh_y**2) # Think about x, y conventions here
-
-    #z = (0. * (r > 1.)) + (21. * (r < 1.))
-    #z = 21. - r
-    #z = (21. - r) * (mesh_x > 0.) * (mesh_y < 0.)
-    z = (21. - r) * np.logical_or(mesh_x > 0., mesh_y > 0.)
-
-    return MaskBand(z, roi)
-    
-############################################################
-
-def readMangleFile(infile, lon, lat, index = None):
-    """
-    DEPRECATED: 2018-05-04
-    Mangle must be set up on your system.
-    The index argument is a temporary file naming convention to avoid file conflicts.
-    Coordinates must be given in the native coordinate system of the Mangle file.
-    """
-    msg = "'mask.readMangleFile': ADW 2018-05-05"
-    DeprecationWarning(msg)
-
-    if index is None:
-        index = np.random.randint(0, 1.e10)
-    
-    coordinate_file = 'temp_coordinate_%010i.dat'%(index)
-    maglim_file = 'temp_maglim_%010i.dat'%(index)
-
-    writer = open(coordinate_file, 'w')
-    for ii in range(0, len(lon)):
-        writer.write('%12.5f%12.5f\n'%(lon[ii], lat[ii]))
-    writer.close()
-
-    os.system('polyid -W %s %s %s || exit'%(infile,
-                                            coordinate_file,
-                                            maglim_file))
-
-    reader = open(maglim_file)
-    lines = reader.readlines()
-    reader.close()
-
-    os.remove(maglim_file)
-    os.remove(coordinate_file)
-
-    maglim = []
-    for ii in range(1, len(lines)):
-        if len(lines[ii].split()) == 3:
-            maglim.append(float(lines[ii].split()[2]))
-        elif len(lines[ii].split()) == 2:
-            maglim.append(0.) # Coordinates outside of the MANGLE ploygon
-        elif len(lines[ii].split()) > 3:
-            msg = 'Coordinate inside multiple polygons, masking that coordinate.'
-            logger.warning(msg)
-            maglim.append(0.)
-        else:
-            msg = 'Cannot parse maglim file, unexpected number of columns.'
-            logger.error(msg)
-            break
-            
-    maglim = np.array(maglim)
-    return maglim
-
-############################################################
-
-def allSkyMask(infile, nside):
-    msg = "'mask.allSkyMask': ADW 2018-05-05"
-    DeprecationWarning(msg)
-    lon, lat = ugali.utils.skymap.allSkyCoordinates(nside)
-    maglim = readMangleFile(infile, lon, lat, index = None)
-    return maglim
-
-############################################################
-
-def scale(mask, mag_scale, outfile=None):
-    """
-    Scale the completeness depth of a mask such that mag_new = mag + mag_scale.
-    Input is a full HEALPix map.
-    Optionally write out the scaled mask as an sparse HEALPix map.
-    """
-    msg = "'mask.scale': ADW 2018-05-05"
-    DeprecationWarning(msg)
-    mask_new = hp.UNSEEN * np.ones(len(mask))
-    mask_new[mask == 0.] = 0.
-    mask_new[mask > 0.] = mask[mask > 0.] + mag_scale
-
-    if outfile is not None:
-        pix = np.nonzero(mask_new > 0.)[0]
-        data_dict = {'MAGLIM': mask_new[pix]}
-        nside = hp.npix2nside(len(mask_new))
-        ugali.utils.skymap.writeSparseHealpixMap(pix, data_dict, nside, outfile)
-
-    return mask_new
-
-############################################################
-
diff --git a/ugali/observation/roi.py b/ugali/observation/roi.py
index adadb25..f79c49f 100644
--- a/ugali/observation/roi.py
+++ b/ugali/observation/roi.py
@@ -18,6 +18,7 @@
 import ugali.utils.projector
 import ugali.utils.skymap
 
+from ugali.utils.logger import logger
 from ugali.utils.config import Config
 from ugali.utils.healpix import query_disc, ang2pix, pix2ang, ang2vec
 
@@ -121,22 +122,12 @@ def __init__(self, config, lon, lat):
         self.delta_mag = self.bins_mag[1] - self.bins_mag[0]
         self.delta_color = self.bins_color[1] - self.bins_color[0]
 
-        ## Axis labels (DEPRECATED)
-        #self.label_x = 'x (deg)'
-        #self.label_y = 'y (deg)'
-        # 
-        #if self.config.params['catalog']['band_1_detection']:
-        #    self.label_mag = '%s (mag)'%(self.config.params['catalog']['mag_1_field'])
-        #else:
-        #    self.label_mag = '%s (mag)'%(self.config.params['catalog']['mag_2_field'])
-        #self.label_color = '%s - %s (mag)'%(self.config.params['catalog']['mag_1_field'],
-        #                                    self.config.params['catalog']['mag_2_field'])
-
     def plot(self, value=None, pixel=None):
         """
         Plot the ROI
         """
-        # DEPRECATED
+        # DEPRECATED: ADW 2021-07-15
+        DeprecationWarning("'roi.plot' is deprecated and will be removed.")
         import ugali.utils.plotting
 
         map_roi = np.array(hp.UNSEEN \
diff --git a/ugali/preprocess/database.py b/ugali/preprocess/database.py
index 13404bf..7894aaf 100755
--- a/ugali/preprocess/database.py
+++ b/ugali/preprocess/database.py
@@ -299,15 +299,15 @@ def footprint(self,nside):
         """
         Download the survey footprint for HEALpix pixels.
         """
-        import healpy
+        import healpy as hp
         import ugali.utils.projector
         if nside > 2**9: raise Exception("Overflow error: nside must be <=2**9")
-        pix = np.arange(healpy.nside2npix(nside),dtype='int')
-        footprint = np.zeros(healpy.nside2npix(nside),dtype='bool')
+        pix = np.arange(hp.nside2npix(nside),dtype='int')
+        footprint = np.zeros(hp.nside2npix(nside),dtype='bool')
         ra,dec = ugali.utils.projector.pixToAng(nside,pix)
         table_name = 'Pix%i'%nside
         self.upload(np.array([pix,ra,dec]).T, ['pix','ra','dec'], name=table_name)
-        radius = healpy.nside2resol(nside_superpix,arcmin=True)
+        radius = hp.nside2resol(nside,arcmin=True)
 
         query="""
         SELECT t.pix, dbo.fInFootprintEq(t.ra, t.dec, %g)
@@ -327,8 +327,7 @@ def __init__(self,release='SVA1_GOLD'):
     def _setup_desdbi(self):
         # Function here to setup trivialAccess client...
         # This should work but it doesn't
-        import warnings
-        warnings.warn("desdbi is deprecated", DeprecationWarning)
+        DeprecationWarning("'desdbi' is deprecated")
         import despydb.desdbi
 
     def generate_query(self, ra_min,ra_max,dec_min,dec_max,filename,db):
diff --git a/ugali/preprocess/padova.py b/ugali/preprocess/padova.py
index 3d4ad03..fd9fa87 100755
--- a/ugali/preprocess/padova.py
+++ b/ugali/preprocess/padova.py
@@ -227,7 +227,7 @@ def verify(cls, filename, survey, age, metallicity):
                 
             try:
                 s = lines[2].split()[-2]
-                assert dict_output[survey][:4] in s
+                assert photname_dict[survey][:4] in s
             except:
                 msg = "Incorrect survey:\n"+lines[2]
                 raise Exception(msg)
diff --git a/ugali/scratch/PlotAllSkyHealpix.py b/ugali/scratch/PlotAllSkyHealpix.py
index 35a871d..fae5ad6 100644
--- a/ugali/scratch/PlotAllSkyHealpix.py
+++ b/ugali/scratch/PlotAllSkyHealpix.py
@@ -1,9 +1,11 @@
 #!/usr/bin/env python
 import healpy
 import pylab as plt
+import numpy
+
 import ugali.utils.skymap
 from ugali.utils.projector import celToGal
-import numpy
+from ugali.utils.logger import logger
 
 default_kwargs = dict( xytext=(5,5),textcoords='offset points',
                        ha="left",va="center",
diff --git a/ugali/scratch/simulation/assemble.py b/ugali/scratch/simulation/assemble.py
index f989fe3..7192b89 100644
--- a/ugali/scratch/simulation/assemble.py
+++ b/ugali/scratch/simulation/assemble.py
@@ -1,3 +1,5 @@
+DeprecationWarning("'assemble.py' should be removed")
+
 import sys
 import os
 import shutil
@@ -9,7 +11,8 @@
 config = yaml.load(open(config_file))
 
 if os.path.exists(outdir):
-    raw_input('Are you sure you want to continue? [Press ENTER to continue]')
+    print("Output directory exists: %s"%outdir)
+    input('Are you sure you want to continue? [Press ENTER to continue]')
 
 if not os.path.exists(outdir):
     os.mkdir(outdir)
diff --git a/ugali/scratch/simulation/merge_population_files.py b/ugali/scratch/simulation/merge_population_files.py
index 25e67c8..1d5553f 100644
--- a/ugali/scratch/simulation/merge_population_files.py
+++ b/ugali/scratch/simulation/merge_population_files.py
@@ -1,9 +1,11 @@
+DeprecationWarning("'merge_population_files.py' should be removed")
+
 import sys
 import glob
 import numpy as np
 import astropy.io.fits as pyfits
 
-print sys.argv
+print(sys.argv)
 
 infiles = sorted(glob.glob(sys.argv[1]))
 outfile = sys.argv[2]
@@ -11,7 +13,7 @@
 data_array = []
 header_array = []
 for infile in infiles:
-    print infile
+    print(infile)
     reader = pyfits.open(infile)
     data_array.append(reader[1].data)
     header_array.append(reader[1].header)
@@ -20,11 +22,11 @@
 data_array = np.concatenate(data_array)
 
 
-print '\nWill write output to %s\n'%(outfile)
+print('\nWill write output to %s\n'%(outfile))
 
 tbhdu = pyfits.BinTableHDU(data_array)
 tbhdu.header = header_array[0]
 
-raw_input('Continue?')
+input('Continue?')
 
 tbhdu.writeto(outfile)
diff --git a/ugali/scratch/simulation/simulate_population.py b/ugali/scratch/simulation/simulate_population.py
index 7cc0847..fd6a364 100755
--- a/ugali/scratch/simulation/simulate_population.py
+++ b/ugali/scratch/simulation/simulate_population.py
@@ -159,7 +159,7 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass,
     flag_too_extended = False
     max_extension = 5.0 # deg
     if a_h >= max_extension:
-        print 'Too extended: a_h = %.2f'%(a_h)
+        print('Too extended: a_h = %.2f'%(a_h))
         a_h = max_extension
         flag_too_extended = True
     # Elliptical kernels take the "extension" as the semi-major axis
@@ -270,7 +270,7 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass,
         pylab.savefig('y3_sat_sim_cmd_%s.png'%('test'), dpi=150.)
         
         print('n_Sigma_p = %i'%(n_sigma_p))
-        raw_input('WAIT')
+        input('WAIT')
         
     satellite=odict(lon=lon[cut_detect], lat=lat[cut_detect], 
                     mag_1=mag_1_meas[cut_detect], mag_2=mag_2_meas[cut_detect], 
@@ -378,8 +378,8 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100,
     mag_extinction_2_array = []
     mc_source_id_array = []
     for ii, mc_source_id in enumerate(population['id']):
-        print 'Simulating satellite (%i/%i) ... mc_source_id = %i'%(ii + 1, n, mc_source_id)
-        print '  distance=%(distance).2e, stellar_mass=%(stellar_mass).2e, r_physical=%(r_physical).2e'%(population[ii])
+        print('Simulating satellite (%i/%i) ... mc_source_id = %i'%(ii + 1, n, mc_source_id))
+        print('  distance=%(distance).2e, stellar_mass=%(stellar_mass).2e, r_physical=%(r_physical).2e'%(population[ii]))
         satellite = catsimSatellite(config, population[ii]['lon'], population[ii]['lat'], 
                                     population[ii]['distance'], population[ii]['stellar_mass'], 
                                     population[ii]['r_physical'],population[ii]['ellipticity'],
@@ -435,7 +435,7 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100,
             mag_extinction_2_array.append(satellite['mag_extinction_2'])
             mc_source_id_array.append(np.tile(mc_source_id, len(satellite['lon'])))
         else:
-            print '  difficulty=%i; satellite not simulated...'%difficulty_population[ii]
+            print('  difficulty=%i; satellite not simulated...'%difficulty_population[ii])
 
     # Concatenate the arrays
     print("Concatenating arrays...")
diff --git a/ugali/scratch/simulation/simulation_match.py b/ugali/scratch/simulation/simulation_match.py
index 866f6c9..7cc8476 100644
--- a/ugali/scratch/simulation/simulation_match.py
+++ b/ugali/scratch/simulation/simulation_match.py
@@ -1,3 +1,5 @@
+DeprecationWarning("'simulation_match.py' should be removed")
+
 import numpy as np
 import astropy.io.fits as pyfits
 import pylab
@@ -40,7 +42,7 @@ def wrap(x):
 data_sim = reader_sim[1].data
 reader_sim.close()
 
-print len(data_sim)
+print(len(data_sim))
 
 """
 match_search, match_sim, angsep = ugali.utils.projector.match(data_search['ra'], data_search['dec'], data_sim['RA'], data_sim['DEC'], tol=1.)
@@ -162,13 +164,13 @@ def wrap(x):
 if save:
     pylab.savefig('sim_match_distance_luminosity.pdf')
 
-print '%20s%20s%20s%20s%20s%20s'%('mc_source_id', 'ra', 'dec', 'distance_modulus', 'fracdet', 'density')
+print('%20s%20s%20s%20s%20s%20s'%('mc_source_id', 'ra', 'dec', 'distance_modulus', 'fracdet', 'density'))
 for index in np.nonzero(cut_why_not)[0]:
-    print '%20i%20.3f%20.3f%20.3f%20.3f%20.3f'%(data_sim['mc_source_id'][index], 
+    print('%20i%20.3f%20.3f%20.3f%20.3f%20.3f'%(data_sim['mc_source_id'][index],
                                                 data_sim['ra'][index], data_sim['dec'][index], 
                                                 data_sim['distance_modulus'][index],
                                                 data_sim['fracdet'][index],
-                                                data_sim['density'][index])
+                                                data_sim['density'][index]))
     
 ############################################################
 
@@ -190,7 +192,7 @@ def wrap(x):
 
 classifier_file = 'trained_classifier.txt'
 if fit:
-    print 'Training the machine learning classifier. This may take a while ...'
+    print('Training the machine learning classifier. This may take a while ...')
     t_start = time.time()
     classifier = sklearn.gaussian_process.GaussianProcessClassifier(1.0 * sklearn.gaussian_process.kernels.RBF(0.5))
     #classifier = sklearn.neighbors.KNeighborsClassifier(3, weights='uniform')
@@ -198,18 +200,18 @@ def wrap(x):
     #classifier = sklearn.svm.SVC(gamma=2, C=1)
     classifier.fit(x[cut_train], y[cut_train])
     t_end = time.time()
-    print '  ... training took %.2f seconds'%(t_end - t_start)
+    print('  ... training took %.2f seconds'%(t_end - t_start))
     # Save the trained classifier
     classifier_data = pickle.dumps(classifier)
     writer = open(classifier_file, 'w')
     writer.write(classifier_data)
     writer.close()
-    print 'Saving machine learning classifier to %s ...'%(classifier_file)
+    print('Saving machine learning classifier to %s ...'%(classifier_file))
     os.system('gzip %s'%(classifier_file))
 else:
-    print 'Loading machine learning classifier from %s ...'%(classifier_file)
+    print('Loading machine learning classifier from %s ...'%(classifier_file))
     if os.path.exists(classifier_file + '.gz') and not os.path.exists(classifier_file):
-        print '  Unzipping...'
+        print('  Unzipping...')
         os.system('gunzip -k %s.gz'%(classifier_file))
     reader = open(classifier_file)
     classifier_data = ''.join(reader.readlines())
diff --git a/ugali/scratch/simulation/survey_selection_function.py b/ugali/scratch/simulation/survey_selection_function.py
index 9719da5..eff06da 100644
--- a/ugali/scratch/simulation/survey_selection_function.py
+++ b/ugali/scratch/simulation/survey_selection_function.py
@@ -1,5 +1,5 @@
 """
-
+Create a survey selection function
 """
 
 import time
@@ -9,6 +9,7 @@
 import numpy as np
 import healpy as hp
 import astropy.io.fits as pyfits
+import matplotlib
 import matplotlib.pyplot as plt
 import pylab
 
@@ -150,7 +151,7 @@ def __init__(self, config_file):
 
     def loadFracdet(self):
         if self.m_fracdet is None:
-            print 'Loading fracdet map from %s ...'%(self.config['infile']['fracdet'])
+            print('Loading fracdet map from %s ...'%(self.config['infile']['fracdet']))
             self.m_fracdet = hp.read_map(self.config['infile']['fracdet'], nest=False)
 
     def loadPopulationMetadata(self):
@@ -165,7 +166,7 @@ def loadSimResults(self):
 
     def loadRealResults(self):
         if self.data_real is None:
-            print 'Loading real data search results from %s ...'%(self.config[self.algorithm]['real_results'])
+            print('Loading real data search results from %s ...'%(self.config[self.algorithm]['real_results']))
             reader = pyfits.open(self.config[self.algorithm]['real_results'])
             self.data_real = reader[1].data
             reader.close()
@@ -206,7 +207,7 @@ def trainClassifier(self):
 
         # Train random forest classifier
         if True:
-            print 'Training the machine learning classifier. This may take a while ...'
+            print('Training the machine learning classifier. This may take a while ...')
             t_start = time.time()
             parameters = {'n_estimators':(500,1000)}#, 'criterion':["gini","entropy"], "min_samples_leaf": [1,2,4]}
             rf = RandomForestClassifier(oob_score=True)
@@ -223,14 +224,14 @@ def trainClassifier(self):
             print(self.classifier.best_estimator_)
             print(self.classifier.best_params_)
             t_end = time.time()
-            print '  ... training took %.2f seconds'%(t_end - t_start)
+            print('  ... training took %.2f seconds'%(t_end - t_start))
 
             # Save the trained classifier
             classifier_data = pickle.dumps(self.classifier)
             writer = open(self.config[self.algorithm]['classifier'], 'w')
             writer.write(classifier_data)
             writer.close()
-            print 'Saving machine learning classifier to %s ...'%(self.config[self.algorithm]['classifier'])
+            print('Saving machine learning classifier to %s ...'%(self.config[self.algorithm]['classifier']))
         else:
             self.loadClassifier()
 
@@ -310,7 +311,8 @@ def validateClassifier(self, cut_detect, cut_train, cut_geometry, y_pred):
                      'actual': 'black',
                      'hsc': 'black'}
 
-        title = 'N_train = %i ; N_test = %i'%(len(cut_train),len(cut_test))        import matplotlib
+        title = 'N_train = %i ; N_detect = %i'%(len(cut_train),len(cut_detect))
+
         cmap = matplotlib.colors.ListedColormap(['Gold', 'Orange', 'DarkOrange', 'OrangeRed', 'Red'])
 
         pylab.figure()
@@ -373,9 +375,9 @@ def validateClassifier(self, cut_detect, cut_train, cut_geometry, y_pred):
         #    pylab.savefig('sim_match_scatter_prediction.pdf')#, dpi=dpi)
 
     def loadClassifier(self):
-        print 'Loading machine learning classifier from %s ...'%(self.config[self.algorithm]['classifier'])
+        print('Loading machine learning classifier from %s ...'%(self.config[self.algorithm]['classifier']))
         if os.path.exists(self.config[self.algorithm]['classifier'] + '.gz') and not os.path.exists(self.config[self.algorithm]['classifier']):
-            print '  Unzipping...'
+            print('  Unzipping...')
             os.system('gunzip -k %s.gz'%(self.config[self.algorithm]['classifier']))
         reader = open(self.config[self.algorithm]['classifier'])
         classifier_data = ''.join(reader.readlines())
@@ -449,7 +451,6 @@ def predict(self, lon, lat, **kwargs):
         return pred, flags_geometry
 
     def validatePredict(self, pred, flags_geometry, lon, lat, r_physical, abs_mag, distance):
-        import matplotlib
         cmap = matplotlib.colors.ListedColormap(['Gold', 'Orange', 'DarkOrange', 'OrangeRed', 'Red'])
 
         pylab.figure()
diff --git a/ugali/scratch/simulation/validate_sim_population.py b/ugali/scratch/simulation/validate_sim_population.py
index 207f384..038ef96 100644
--- a/ugali/scratch/simulation/validate_sim_population.py
+++ b/ugali/scratch/simulation/validate_sim_population.py
@@ -48,11 +48,11 @@ def getCatalog(catalog_dir):
     catalog_infiles = sorted(glob.glob(catalog_dir + '/*catalog*.fits'))
     data_array = []
     for catalog_infile in catalog_infiles:
-        print '  Reading %s ...'%(catalog_infile)
+        print('  Reading %s ...'%(catalog_infile))
         reader = pyfits.open(catalog_infile)
         data_array.append(reader[1].data)
         reader.close()
-    print '  Merging ...'
+    print('  Merging ...')
     return np.concatenate(data_array)
 
 ##########
@@ -64,7 +64,7 @@ def getCatalog(catalog_dir):
 infile_population = 'v4/sim_population_v4.fits'
 reader_population = pyfits.open(infile_population)
 data_population = reader_population[1].data
-print len(data_population)
+print(len(data_population))
 data_population = data_population #[0:500]
 reader_population.close()
 
@@ -207,7 +207,7 @@ def getCatalog(catalog_dir):
 """
 ##########
 """
-print "Machine learning"
+print("Machine learning")
 
 save = False
 
diff --git a/ugali/simulation/analyzer.py b/ugali/simulation/analyzer.py
index 1d0f8d3..a18b4d0 100755
--- a/ugali/simulation/analyzer.py
+++ b/ugali/simulation/analyzer.py
@@ -30,16 +30,13 @@
 from ugali.utils import mlab
 
 # Analysis flags
-FLAGS = odict([
-        ('FLAG_PROC'    , 0 ),    # Simulation was processed
-        ('FLAG_NOPROC'  , 1 ),    # No processing
-        ('FLAG_NOBJ'    , 2 ),    # Too many catalog objects
-        ('FLAG_FIT'     , 4 ),    # Fit failure
-        ('FLAG_EBV'     , 8 ),    # EBV value too large
-        ('FLAG_MEM'     , 16),    # Memory error
-        ])
-for k,v in FLAGS.items():
-    globals()[k] = v
+FLAGS = odict([])
+FLAGS['FLAG_PROC'  ] = FLAG_PROC   = 0  # Simulation was processed
+FLAGS['FLAG_NOPROC'] = FLAG_NOPROC = 1  # No processing
+FLAGS['FLAG_NOBJ'  ] = FLAG_NOBJ   = 2  # Too many catalog objects
+FLAGS['FLAG_FIT'   ] = FLAG_FIT    = 4  # Fit failure
+FLAGS['FLAG_EBV'   ] = FLAG_EBV    = 8  # EBV value too large
+FLAGS['FLAG_MEM'   ] = FLAG_MEM    = 16 # Memory error
 
 def update_header_flags(filename):
     fits = fitsio.FITS(filename,'rw')
diff --git a/ugali/simulation/population.py b/ugali/simulation/population.py
index 0bc4bf7..1e8a9b4 100644
--- a/ugali/simulation/population.py
+++ b/ugali/simulation/population.py
@@ -104,95 +104,6 @@ def satellitePopulation(mask, nside_pix, size,
 
 ############################################################
 
-# ADW: 2019-09-01 DEPRECATED
-def satellitePopulationOrig(config, n,
-                            range_distance_modulus=[16.5, 24.],
-                            range_stellar_mass=[1.e2, 1.e5],
-                            range_r_physical=[5.e-3, 1.],
-                            mode='mask',
-                            plot=False):
-    """
-    Create a population of `n` randomly placed satellites within a
-    survey mask or catalog specified in the config file.  Satellites
-    are distributed uniformly in distance modulus, uniformly in
-    log(stellar_mass / Msun), and uniformly in
-    log(r_physical/kpc). The ranges can be set by the user.
-
-    Returns the simulated area (deg^2) as well as the
-    Parameters
-    ----------
-    config : configuration file or object
-    n      : number of satellites to simulate
-    range_distance_modulus : range of distance modulus to sample
-    range_stellar_mass : range of stellar mass to sample (Msun)
-    range_r_physical : range of azimuthally averaged half light radius (kpc)
-    mode : how to sample the area['mask','catalog']
-
-    Returns
-    -------
-    lon (deg), lat (deg), distance modulus, stellar mass (Msun), and
-    half-light radius (kpc) for each satellite
-    """
-    
-    if type(config) == str:
-        config = ugali.utils.config.Config(config)
-
-    if mode == 'mask':
-        mask_1 = ugali.utils.skymap.readSparseHealpixMap(config.params['mask']['infile_1'], 'MAGLIM')
-        mask_2 = ugali.utils.skymap.readSparseHealpixMap(config.params['mask']['infile_2'], 'MAGLIM')
-        input = (mask_1 > 0.) * (mask_2 > 0.)
-    elif mode == 'catalog':
-        catalog = ugali.observation.catalog.Catalog(config)
-        input = np.array([catalog.lon, catalog.lat])
-    
-    lon, lat, simulation_area = ugali.utils.skymap.randomPositions(input,
-                                                                   config.params['coords']['nside_likelihood_segmentation'],
-                                                                   n=n)
-    distance_modulus = np.random.uniform(range_distance_modulus[0], 
-                                         range_distance_modulus[1], 
-                                         n)
-    stellar_mass = 10**np.random.uniform(np.log10(range_stellar_mass[0]), 
-                                         np.log10(range_stellar_mass[1]), 
-                                         n)
-    
-    half_light_radius_physical = 10**np.random.uniform(np.log10(range_half_light_radius_physical[0]), 
-                                                       np.log10(range_half_light_radius_physical[0]), 
-                                                       n) # kpc
-
-    half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \
-                                             / ugali.utils.projector.distanceModulusToDistance(distance_modulus)))
-    
-    # One choice of theory prior
-    #half_light_radius_physical = ugali.analysis.kernel.halfLightRadius(stellar_mass) # kpc
-    #half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \
-    #                                               / ugali.utils.projector.distanceModulusToDistance(distance_modulus)))
-
-    if plot:
-        pylab.figure()
-        #pylab.scatter(lon, lat, c=distance_modulus, s=500 * half_light_radius)
-        #pylab.colorbar()
-        pylab.scatter(lon, lat, edgecolors='none')
-        xmin, xmax = pylab.xlim() # Reverse azimuthal axis
-        pylab.xlim([xmax, xmin])
-        pylab.title('Random Positions in Survey Footprint')
-        pylab.xlabel('Longitude (deg)')
-        pylab.ylabel('Latitude (deg)')
-
-        pylab.figure()
-        pylab.scatter(stellar_mass, ugali.utils.projector.distanceModulusToDistance(distance_modulus),
-                      c=(60. * half_light_radius), s=500 * half_light_radius, edgecolors='none')
-        pylab.xscale('log')
-        pylab.yscale('log')
-        pylab.xlim([0.5 * range_stellar_mass[0], 2. * range_stellar_mass[1]])
-        pylab.colorbar()
-        pylab.title('Half-light Radius (arcmin)')
-        pylab.xlabel('Stellar Mass (arcmin)')
-        pylab.ylabel('Distance (kpc)')
-
-    return simulation_area, lon, lat, distance_modulus, stellar_mass, half_light_radius 
-
-############################################################
-
 def interpolate_absolute_magnitude():
     iso = ugali.isochrone.factory('Bressan2012',age=12,z=0.00010)
 
@@ -272,26 +183,3 @@ def knownPopulation(dwarfs, mask, nside_pix, size):
     population = np.hstack(results)
     population['id'] = np.arange(size)
     return area, population
-
-def plot_population(population):
-    # ADW: DEPRECATED: 2019-09-01
-    pylab.figure()
-    #pylab.scatter(lon, lat, c=distance_modulus, s=500 * half_light_radius)
-    #pylab.colorbar()
-    pylab.scatter(lon, lat, edgecolors='none')
-    xmin, xmax = pylab.xlim() # Reverse azimuthal axis
-    pylab.xlim([xmax, xmin])
-    pylab.title('Random Positions in Survey Footprint')
-    pylab.xlabel('Longitude (deg)')
-    pylab.ylabel('Latitude (deg)')
-
-    pylab.figure()
-    pylab.scatter(stellar_mass, distance,c=r_physical,
-                  s=500 * r_physical, edgecolors='none')
-    pylab.xscale('log')
-    pylab.yscale('log')
-    pylab.xlim([0.5 * range_stellar_mass[0], 2. * range_stellar_mass[1]])
-    pylab.colorbar()
-    pylab.title('Half-light Radius (arcmin)')
-    pylab.xlabel('Stellar Mass (arcmin)')
-    pylab.ylabel('Distance (kpc)')
diff --git a/ugali/utils/batch.py b/ugali/utils/batch.py
index e8cb4d7..371c240 100644
--- a/ugali/utils/batch.py
+++ b/ugali/utils/batch.py
@@ -312,7 +312,7 @@ def parse_options(self, **opts):
 
 class Condor(Batch):
     """ Not implemented yet... """
-    def __init__(self):
+    def __init__(self, **kwargs):
         super(Condor,self).__init__(**kwargs)
         logger.warning('Condor cluster is untested')
         
diff --git a/ugali/utils/config.py b/ugali/utils/config.py
index 6a0d17e..535dac2 100644
--- a/ugali/utils/config.py
+++ b/ugali/utils/config.py
@@ -200,9 +200,6 @@ def _createFilenames(self,pixels=None):
         """
         nside_catalog = self['coords']['nside_catalog']
 
-        # Deprecated: ADW 2018-06-17
-        #if nside_catalog is None:
-        #    pixels = [None]
         if pixels is not None:
             pixels = [pixels] if np.isscalar(pixels) else pixels
         else:
@@ -225,6 +222,7 @@ def _createFilenames(self,pixels=None):
             if pix is None:
                 # DEPRECTATED: ADW 2018-06-17
                 # This is not really being used anymore
+                raise ValueError('pix cannot be None')
                 catalog = os.path.join(catalog_dir,catalog_base)
                 mask_1 = os.path.join(mask_dir,mask_base_1)
                 mask_2 = os.path.join(mask_dir,mask_base_2)
diff --git a/ugali/utils/fileio.py b/ugali/utils/fileio.py
index 9f1caa1..04708fa 100644
--- a/ugali/utils/fileio.py
+++ b/ugali/utils/fileio.py
@@ -66,7 +66,23 @@ def write(filename,data,**kwargs):
 
     msg = "Unrecognized file type: %s"%filename
     raise ValueError(msg)
-    
+
+def check_formula(formula):
+    """ Check that a formula is valid. """
+    if not (('data[' in formula) and (']' in formula)):
+        msg = 'Invalid formula:\n %s'%formula
+        raise ValueError(msg)
+
+def parse_formula(formula):
+    """ Return the columns used in a formula. """
+    check_formula(formula)
+
+    columns = []
+    for x in formula.split('data[')[1:]:
+        col = x.split(']')[0].replace('"','').replace("'",'')
+        columns += [col]
+
+    return columns
 
 def add_column(filename,column,formula,force=False):
     """ Add a column to a FITS file.
@@ -246,65 +262,6 @@ def insert_columns(filename,data,ext=1,force=False,colnum=None):
 # Dealing with FITS files
 def write_fits(filename,data,header=None,force=False):
     if os.path.exists(filename) and not force:
-        found(filename)
+        logger.found(filename)
         return
     fitsio.write(filename,data,header=header,clobber=force)
-
-# Writing membership files
-def write_membership(loglike,filename):
-    """
-    Write a catalog file of the likelihood region including
-    membership properties.
-
-    Parameters:
-    -----------
-    loglike : input loglikelihood object
-    filename : output filename
-    
-    Returns:
-    --------
-    None
-    """
-
-    ra,dec = gal2cel(loglike.catalog.lon,loglike.catalog.lat)
-        
-    name_objid = loglike.config['catalog']['objid_field']
-    name_mag_1 = loglike.config['catalog']['mag_1_field']
-    name_mag_2 = loglike.config['catalog']['mag_2_field']
-    name_mag_err_1 = loglike.config['catalog']['mag_err_1_field']
-    name_mag_err_2 = loglike.config['catalog']['mag_err_2_field']
-
-    # Angular and isochrone separations
-    sep = angsep(loglike.source.lon,loglike.source.lat,
-                 loglike.catalog.lon,loglike.catalog.lat)
-    isosep = loglike.isochrone.separation(loglike.catalog.mag_1,loglike.catalog.mag_2)
-
-    data = odict()
-    data[name_objid] = loglike.catalog.objid
-    data['GLON'] = loglike.catalog.lon
-    data['GLAT'] = loglike.catalog.lat
-    data['RA']   = ra
-    data['DEC']  = dec
-    data[name_mag_1] = loglike.catalog.mag_1
-    data[name_mag_err_1] = loglike.catalog.mag_err_1
-    data[name_mag_2] = loglike.catalog.mag_2
-    data[name_mag_err_2] = loglike.catalog.mag_err_2
-    data['COLOR'] = loglike.catalog.color
-    data['ANGSEP'] = sep
-    data['ISOSEP'] = isosep
-    data['PROB'] = loglike.p
-
-    # HIERARCH allows header keywords longer than 8 characters
-    header = []
-    for param,value in loglike.source.params.items():
-        card = dict(name='HIERARCH %s'%param.upper(),
-                    value=value.value,
-                    comment=param)
-        header.append(card)
-    card = dict(name='HIERARCH %s'%'TS',value=loglike.ts(),
-                comment='test statistic')
-    header.append(card)
-    card = dict(name='HIERARCH %s'%'TIMESTAMP',value=time.asctime(),
-                comment='creation time')
-    header.append(card)
-    fitsio.write(filename,data,header=header,clobber=True)
diff --git a/ugali/utils/healpix.py b/ugali/utils/healpix.py
index d92b4aa..e4e721b 100644
--- a/ugali/utils/healpix.py
+++ b/ugali/utils/healpix.py
@@ -8,11 +8,12 @@
 import numpy
 import numpy as np
 import healpy as hp
-import healpy
+import fitsio
 
 import ugali.utils.projector
 from ugali.utils import fileio
 from ugali.utils.logger import logger
+import ugali.utils.fileio
 
 ############################################################
 
@@ -382,7 +383,6 @@ def write_partial_map(filename, data, nside, coord=None, nest=False,
     --------
     None
     """
-    import fitsio
 
     # ADW: Do we want to make everything uppercase?
 
@@ -590,15 +590,15 @@ def read_map(filename, nest=False, hdu=None, h=False, verbose=True):
     m [, header] : array, optionally with header appended
       The map read from the file, and the header if *h* is True.
     """
-    
+    import fitsio
     data,hdr = fitsio.read(filename,header=True,ext=hdu)
 
     nside = int(hdr.get('NSIDE'))
     if verbose: print('NSIDE = {0:d}'.format(nside))
 
-    if not healpy.isnsideok(nside):
+    if not hp.isnsideok(nside):
         raise ValueError('Wrong nside parameter.')
-    sz=healpy.nside2npix(nside)
+    sz=hp.nside2npix(nside)
 
     ordering = hdr.get('ORDERING','UNDEF').strip()
     if verbose: print('ORDERING = {0:s} in fits file'.format(ordering))
@@ -615,28 +615,27 @@ def read_map(filename, nest=False, hdu=None, h=False, verbose=True):
 
     # Could be done more efficiently (but complicated) by reordering first
     if hdr['INDXSCHM'] == 'EXPLICIT':
-        m = healpy.UNSEEN*np.ones(sz,dtype=data[fields[1]].dtype)
+        m = hp.UNSEEN*np.ones(sz,dtype=data[fields[1]].dtype)
         m[data[fields[0]]] = data[fields[1]]
     else:
         m = data[fields[0]].ravel()
 
-    if (not healpy.isnpixok(m.size) or (sz>0 and sz != m.size)) and verbose:
+    if (not hp.isnpixok(m.size) or (sz>0 and sz != m.size)) and verbose:
         print('nside={0:d}, sz={1:d}, m.size={2:d}'.format(nside,sz,m.size))
         raise ValueError('Wrong nside parameter.')
     if nest is not None:
         if nest and ordering.startswith('RING'):
-            idx = healpy.nest2ring(nside,np.arange(m.size,dtype=np.int32))
+            idx = hp.nest2ring(nside,np.arange(m.size,dtype=np.int32))
             m = m[idx]
             if verbose: print('Ordering converted to NEST')
         elif (not nest) and ordering.startswith('NESTED'):
-            idx = healpy.ring2nest(nside,np.arange(m.size,dtype=np.int32))
+            idx = hp.ring2nest(nside,np.arange(m.size,dtype=np.int32))
             m = m[idx]
             if verbose: print('Ordering converted to RING')
 
     if h: return m, hdr
     else: return m
 
-
 if __name__ == "__main__":
     import argparse
     description = __doc__
diff --git a/ugali/utils/idl.py b/ugali/utils/idl.py
index d6efdcf..14bced8 100644
--- a/ugali/utils/idl.py
+++ b/ugali/utils/idl.py
@@ -97,8 +97,8 @@ def jprecess(ra, dec, mu_radec=None, parallax=None, rad_vel=None, epoch=None):
       ra = array(ra)
       dec = array(dec)
    else:
-      ra = array([ra0])
-      dec = array([dec0])
+      ra = array([ra])
+      dec = array([dec])
    n = ra.size
     
    if rad_vel is None:   
@@ -111,7 +111,7 @@ def jprecess(ra, dec, mu_radec=None, parallax=None, rad_vel=None, epoch=None):
     
    if (mu_radec is not None):   
       if (array(mu_radec).size != 2 * n):   
-         raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,' + strtrim(n, 2) + ')')
+         raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,%s)'%n)
       mu_radec = mu_radec * 1.
     
    if parallax is None:   
@@ -315,7 +315,7 @@ def bprecess(ra0, dec0, mu_radec=None, parallax=None, rad_vel=None, epoch=None):
    
    if (mu_radec is not None):   
       if (array(mu_radec).size != 2 * n):   
-         raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,' + strtrim(n, 2) + ')')
+         raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,%s)'%n)
       mu_radec = mu_radec * 1.
    
    if parallax is None:   
diff --git a/ugali/utils/logger.py b/ugali/utils/logger.py
index 8a2d454..8863abc 100644
--- a/ugali/utils/logger.py
+++ b/ugali/utils/logger.py
@@ -3,6 +3,7 @@
 Interface to python logging. For more info see:
 http://docs.python.org/2/howto/logging.html
 """
+import os.path
 import logging
 
 class SpecialFormatter(logging.Formatter):
diff --git a/ugali/utils/plotting.py b/ugali/utils/plotting.py
index ff2953b..c289bf3 100644
--- a/ugali/utils/plotting.py
+++ b/ugali/utils/plotting.py
@@ -585,8 +585,6 @@ def drawCMD(self, ax=None, radius=None, zidx=None):
         ax.set_xlabel('Color (mag)')
         ax.set_ylabel('Magnitude (mag)')
 
-        try:    ax.cax.colorbar(im)
-        except: pass
         ax.annotate("Stars",**self.label_kwargs)
 
 
@@ -1417,15 +1415,6 @@ def plot_chain(chain,burn=None,clip=None):
 
 ###################################################
 
-    
-def plotSkymapCatalog(lon,lat,**kwargs):
-    """
-    Plot a catalog of coordinates on a full-sky map.
-    """
-    fig = plt.figure()
-    ax = plt.subplot(111,projection=projection)
-    drawSkymapCatalog(ax,lon,lat,**kwargs)
-
 def drawSkymapCatalog(ax,lon,lat,**kwargs):
     mapping = {
         'ait':'aitoff',
diff --git a/ugali/utils/projector.py b/ugali/utils/projector.py
index 4e6f838..75e9b7b 100644
--- a/ugali/utils/projector.py
+++ b/ugali/utils/projector.py
@@ -4,7 +4,7 @@
 Based on Calabretta & Greisen 2002, A&A, 357, 1077-1122
 http://adsabs.harvard.edu/abs/2002A%26A...395.1077C
 """
-
+import re
 import numpy as np
 
 from ugali.utils.logger import logger
@@ -31,8 +31,9 @@ def setReference(self, lon_ref, lat_ref, zenithal=False):
         if not zenithal:
             phi = (-np.pi / 2.) + np.radians(lon_ref)
             theta = np.radians(lat_ref)
-            psi = np.radians(90.) # psi = 90 corresponds to (0, 0) psi = -90 corresponds to (180, 0)
-        
+            # psi = 90 corresponds to (0, 0)
+            # psi = -90 corresponds to (180, 0)
+            psi = np.radians(90.)
 
         cos_psi,sin_psi = np.cos(psi),np.sin(psi)
         cos_phi,sin_phi = np.cos(phi),np.sin(phi)
@@ -449,9 +450,10 @@ def hms2dec(hms):
     return decimal
 
 def dms2dec(dms):
-    """
-    Convert latitude from degrees,minutes,seconds in string or 3-array
-   format to decimal degrees.
+    """Convert latitude from degrees,minutes,seconds in string or 3-array
+    format to decimal degrees.
+
+    ADW: This really should be replaced by astropy
     """
     DEGREE = 360.
     HOUR = 24.
@@ -463,7 +465,7 @@ def dms2dec(dms):
     # http://docs.scipy.org/doc/numpy-1.7.0/reference/c-api.coremath.html#NPY_NZERO
 
     if isstring(dms):
-        degree,minute,second = np.array(re.split('[dms]',hms))[:3].astype(float)
+        degree,minute,second = np.array(re.split('[dms]',dms))[:3].astype(float)
     else:
         degree,minute,second = dms.T
 
diff --git a/ugali/utils/skymap.py b/ugali/utils/skymap.py
index 5937ff3..cca24e8 100644
--- a/ugali/utils/skymap.py
+++ b/ugali/utils/skymap.py
@@ -84,23 +84,6 @@ def inFootprint(config, pixels, nside=None):
     
     return inside
 
-def footprint(config, nside=None):
-    """
-    UNTESTED.
-    Should return a boolean array representing the pixels in the footprint.
-    """
-    DeprecationWarning("This function is deprecated.")
-    config = Config(config)
-    if nside is None:
-        nside = config['coords']['nside_pixel']
-    elif nside < config['coords']['nside_catalog']:
-        raise Exception('Requested nside=%i is greater than catalog_nside'%nside)
-    elif nside > config['coords']['nside_pixel']:
-        raise Exception('Requested nside=%i is less than pixel_nside'%nside)
-    pix = np.arange(hp.nside2npix(nside), dtype=int)
-    return inFootprint(config,pix)
-
-
 ############################################################
 
 def allSkyCoordinates(nside):
diff --git a/ugali/utils/stats.py b/ugali/utils/stats.py
index 7e41e90..aec8a15 100644
--- a/ugali/utils/stats.py
+++ b/ugali/utils/stats.py
@@ -4,6 +4,7 @@
 """
 
 import copy
+from collections import OrderedDict as odict
 
 import numpy as np
 import numpy.lib.recfunctions as recfuncs