diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index fa862a1..90d992b 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:
@@ -32,25 +32,21 @@ jobs:
         conda info
     - name: Create conda environment
       run: |
-        conda create -q -n env python=${{ matrix.python-version }} numpy scipy matplotlib astropy healpy pyyaml emcee fitsio corner -c conda-forge -c kadrlica
+        conda create -y -q -n env python=${{ matrix.python-version }} numpy scipy matplotlib astropy healpy pyyaml emcee fitsio corner -c conda-forge -c kadrlica
+        # Add UGALIDIR to environment
+        conda env config vars set UGALIDIR=$HOME/.ugali -n env
     - name: Install package
       run: |
         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
+        conda env config vars list
+        python setup.py -q install --isochrones --catalogs --tests
     - 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: Test with nose
       run: |
         source activate env
diff --git a/README.md b/README.md
index 2b1adc5..740f056 100644
--- a/README.md
+++ b/README.md
@@ -50,7 +50,7 @@ The `ugali` source code is distributed with several auxiliary libraries for isoc
 
 ```
 cd $UGALIDIR
-wget https://github.com/kadrlica/ugali/releases/download/v1.7.0rc0/ugali-des-bressan2012.tar.gz
+wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.8.0/ugali-des-bressan2012.tar.gz
 tar -xzf ugali-des-bressan2012.tar.gz
 ```
 
diff --git a/setup.py b/setup.py
index 201c725..217cf1f 100644
--- a/setup.py
+++ b/setup.py
@@ -24,25 +24,27 @@ def find_packages():
 HERE = os.path.abspath(os.path.dirname(__file__))
 URL = 'https://github.com/DarkEnergySurvey/ugali'
 DESC = "Ultra-faint galaxy likelihood toolkit."
-LONG_DESC = "See %s"%URL
+LONG_DESC = "%s\n%s"%(DESC,URL)
 CLASSIFIERS = """\
-Development Status :: 3 - Alpha
+Development Status :: 4 - Beta
 Intended Audience :: Science/Research
 Intended Audience :: Developers
 License :: OSI Approved :: MIT License
 Natural Language :: English
 Operating System :: MacOS :: MacOS X
 Operating System :: POSIX :: Linux
-Programming Language :: Python :: 2.7
+Programming Language :: Python :: 2
+Programming Language :: Python :: 3
 Topic :: Scientific/Engineering
 Topic :: Scientific/Engineering :: Astronomy
 Topic :: Scientific/Engineering :: Physics
 """
 
-RELEASE = URL+'/releases/download/v1.7.0'
+RELEASE_URL = URL+'/releases/download/v1.8.0'
 UGALIDIR = os.getenv("UGALIDIR","$HOME/.ugali")
 ISOSIZE = "~1MB" 
 CATSIZE = "~20MB"
+TSTSIZE = "~1MB"
 # Could find file size dynamically, but it's a bit slow...
 # int(urllib.urlopen(ISOCHRONES).info().getheaders("Content-Length")[0])/1024**2
 SURVEYS = ['des','ps1','sdss','lsst']
@@ -62,7 +64,7 @@ def read(self, size):
     def progress_bar(count, block_size, total_size):
         block = 100*block_size/float(total_size)
         progress = count*block
-        if progress % 1 < 1.01*block:
+        if progress % 5 < 1.01*block:
             msg = '\r[{:51}] ({:d}%)'.format(int(progress//2)*'='+'>',int(progress))
             sys.stdout.write(msg)
             sys.stdout.flush()
@@ -77,7 +79,7 @@ class TarballCommand(distutils.cmd.Command,object):
          'force installation (overwrite any existing files)')
         ]
     boolean_options = ['force']
-    release = RELEASE
+    release_url = RELEASE_URL
     _tarball = None
     _dirname = None
 
@@ -114,7 +116,7 @@ def install_tarball(self, tarball):
             os.makedirs(self.ugali_dir)
         os.chdir(self.ugali_dir)
 
-        url = os.path.join(self.release,tarball)
+        url = os.path.join(self.release_url,tarball)
 
         print("downloading %s..."%url)
         if urlopen(url).getcode() >= 400:
@@ -158,6 +160,12 @@ class CatalogCommand(TarballCommand):
     _tarball = 'ugali-catalogs.tar.gz'
     _dirname = 'catalogs'
 
+class TestsCommand(TarballCommand):
+    """ Command for downloading catalog files """
+    description = "install test data"
+    _tarball = 'ugali-test-data.tar.gz'
+    _dirname = 'testdata'
+
 class IsochroneCommand(TarballCommand):
     """ Command for downloading isochrone files """
     description = "install isochrone files"
@@ -223,6 +231,7 @@ class install(_install):
     user_options = _install.user_options + [
         ('isochrones',None,"install isochrone files (%s)"%ISOSIZE),
         ('catalogs',None,"install catalog files (%s)"%CATSIZE),
+        ('tests',None,"install test data (%s)"%TSTSIZE),
         ('ugali-dir=',None,"install file directory [default: %s]"%UGALIDIR),
     ]
     boolean_options = _install.boolean_options + ['isochrones','catalogs']
@@ -232,6 +241,7 @@ def initialize_options(self):
         self.ugali_dir = os.path.expandvars(UGALIDIR)
         self.isochrones = False
         self.catalogs = False
+        self.tests = False
 
     def run(self):
         # run superclass install
@@ -246,7 +256,10 @@ def run(self):
 
         if self.catalogs: 
             self.install_catalogs()
-            
+
+        if self.tests:
+            self.install_tests()
+
     def install_isochrones(self):
         """
         Call to isochrone install command:
@@ -266,10 +279,21 @@ def install_catalogs(self):
         cmd_obj.force = self.force
         if self.ugali_dir: cmd_obj.ugali_dir = self.ugali_dir
         self.run_command('catalogs')
-            
+
+    def install_tests(self):
+        """
+        Call to catalog install command:
+        http://stackoverflow.com/a/24353921/4075339
+        """
+        cmd_obj = self.distribution.get_command_obj('tests')
+        cmd_obj.force = self.force
+        if self.ugali_dir: cmd_obj.ugali_dir = self.ugali_dir
+        self.run_command('tests')
+
 CMDCLASS = versioneer.get_cmdclass()
 CMDCLASS['isochrones'] = IsochroneCommand
 CMDCLASS['catalogs'] = CatalogCommand
+CMDCLASS['tests'] = TestsCommand
 CMDCLASS['install'] = install
 
 setup(
@@ -278,9 +302,11 @@ def install_catalogs(self):
     cmdclass=CMDCLASS,
     url=URL,
     author='Keith Bechtol & Alex Drlica-Wagner',
-    author_email='bechtol@kicp.uchicago.edu, kadrlica@fnal.gov',
+    author_email='bechtol@wisc.edu, kadrlica@fnal.gov',
     scripts = [],
     install_requires=[
+        'astropy',
+        'matplotlib',
         'numpy >= 1.9.0',
         'scipy >= 0.14.0',
         'healpy >= 1.6.0',
@@ -288,7 +314,6 @@ def install_catalogs(self):
         'emcee >= 2.1.0',
         'corner >= 1.0.0',
         'pyyaml >= 3.10',
-        # Add astropy, fitsio, matplotlib, ...
     ],
     packages=find_packages(),
     description=DESC,
diff --git a/tests/config.yaml b/tests/config.yaml
index b2278e3..e8116e0 100644
--- a/tests/config.yaml
+++ b/tests/config.yaml
@@ -27,7 +27,7 @@ coords:
 
 catalog:
   # Color always defined as mag_1 - mag_2
-  dirname: ./healpix                    # directory of catalog files
+  dirname: ${UGALIDIR}/testdata/healpix          # directory of catalog files
   basename: "catalog_hpx%04i.fits"      # catalog file basename format
   objid_field       : COADD_OBJECT_ID   # unique object identifier
   lon_field         : RA                # name of longitude field
@@ -45,7 +45,7 @@ catalog:
   selection         : "(np.abs(self.data['WAVG_SPREAD_MODEL_I']) < 0.003)"
 
 data: 
-  dirname: ./raw
+  dirname: ${UGALIDIR}/raw
   script : ./ugali/preprocess/database.py
   survey : des
   release: y3a2
@@ -53,16 +53,11 @@ data:
   footprint: ./maps/footprint_nside4096_equ.fits.gz
 
 mask:
-  dirname    : ./mask
+  dirname    : ${UGALIDIR}/testdata/mask
   basename_1 : "maglim_g_hpx%04i.fits"
   basename_2 : "maglim_r_hpx%04i.fits"
   minimum_solid_angle: 0.1 # deg^2
 
-mangle:
-  dirname    : ./maps
-  filename_1 : 'y3a2_g_o.4096_t.32768_maglim_EQU.fits.gz'
-  filename_2 : 'y3a2_r_o.4096_t.32768_maglim_EQU.fits.gz'
-
 #ADW: Depricated in favor of 'binning'
 color:
   min   : &cmin -0.5
diff --git a/tests/test_config.py b/tests/test_config.py
index 172b993..305ea28 100644
--- a/tests/test_config.py
+++ b/tests/test_config.py
@@ -15,12 +15,12 @@ def test_config():
     np.testing.assert_equal(len(config.filenames),768)
 
     np.testing.assert_equal(config.filenames['pix'].compressed()[0],687)
-    np.testing.assert_equal(config.filenames['catalog'].compressed()[0],
-                            './healpix/catalog_hpx0687.fits')
-    np.testing.assert_equal(config.filenames['mask_1'].compressed()[0],
-                            './mask/maglim_g_hpx0687.fits')
-    np.testing.assert_equal(config.filenames['mask_2'].compressed()[0],
-                            './mask/maglim_r_hpx0687.fits')
+    catfile = os.path.basename(config.filenames['catalog'].compressed()[0])
+    np.testing.assert_equal(catfile,'catalog_hpx0687.fits')
+    maskfile = os.path.basename(config.filenames['mask_1'].compressed()[0])
+    np.testing.assert_equal(maskfile,'maglim_g_hpx0687.fits')
+    maskfile = os.path.basename(config.filenames['mask_2'].compressed()[0])
+    np.testing.assert_equal(maskfile,'maglim_r_hpx0687.fits')
 
     np.testing.assert_equal(config.likefile,'./scan/scan_%08i_%s.fits')
     np.testing.assert_equal(config.mergefile,'./scan/merged_scan.fits')
diff --git a/tests/test_observation.py b/tests/test_observation.py
index e97887b..c9e136a 100644
--- a/tests/test_observation.py
+++ b/tests/test_observation.py
@@ -44,8 +44,13 @@ def test_catalog():
     """ Test ugali.observation.catalog """
     import ugali.observation.roi
     import ugali.observation.catalog
+    import ugali.utils.config
 
-    filename='healpix/catalog_hpx0687.fits'
+    # Get catalog filename
+    config = ugali.utils.config.Config(CONFIG)
+    filename = config.filenames['catalog'].compressed()[0]
+
+    # Create catalog from filename
     catalog = ugali.observation.catalog.Catalog(CONFIG,filenames=filename)
 
     roi = ugali.observation.roi.ROI(CONFIG, LON, LAT)
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 1e8fc36..944cbf9 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..0b05fed 100644
--- a/ugali/isochrone/mesa.py
+++ b/ugali/isochrone/mesa.py
@@ -79,7 +79,7 @@ class Dotter2016(Isochrone):
             des = odict([
                 (2, ('mass_init',float)),
                 (3, ('mass_act',float)),
-                (8, ('log_lum',float)),
+                (6, ('log_lum',float)),
                 (9, ('u',float)),
                 (10,('g',float)),
                 (11,('r',float)),
@@ -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 2ac1af9..d673934 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 7890ff4..77b03cb 100644
--- a/ugali/isochrone/parsec.py
+++ b/ugali/isochrone/parsec.py
@@ -418,7 +418,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...
@@ -503,6 +503,7 @@ class Marigo2017(ParsecIsochrone):
                 (30,('Y',float)),
                 ]),
         )
+
     columns['lsst'] = copy.deepcopy(columns['lsst_dp0'])
 
     def _find_column_numbers(self):
@@ -526,7 +527,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..611f5ff 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
@@ -627,7 +592,7 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None):
         index_color = np.arange(len(self.roi.centers_color))
         # Add the cumulative leakage back into the last bin of the CMD
         leakage = (cmd_background * ~observable).sum(axis=0)
-        cmd_background[[index_mag,index_color]] += leakage
+        cmd_background[(index_mag,index_color)] += leakage
         # Zero out all non-observable bins
         cmd_background *= observable
 
@@ -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/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/analysis/color_lut.py b/ugali/scratch/color_lut.py
similarity index 100%
rename from ugali/analysis/color_lut.py
rename to ugali/scratch/color_lut.py
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 d88286a..a64623b 100755
--- a/ugali/scratch/simulation/simulate_population.py
+++ b/ugali/scratch/simulation/simulate_population.py
@@ -290,7 +290,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], 
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..78d62de 100644
--- a/ugali/utils/config.py
+++ b/ugali/utils/config.py
@@ -172,88 +172,6 @@ def write(self, filename):
             raise Exception('Unrecognized config format: %s'%ext)
         writer.close()
 
-    def _createFilenames(self,pixels=None):
-        """
-        Create a masked records array of all filenames for the given set of
-        pixels and store the existence of those files in the mask values.
-
-        Examples:
-        f = getFilenames([1,2,3])
-        # All possible catalog files
-        f['catalog'].data
-        # All existing catalog files
-        f['catalog'][~f.mask['catalog']]
-        # or
-        f['catalog'].compressed()
-        # All missing mask_1 files
-        f['mask_1'][f.mask['mask_1']]
-        # Pixels where all files exist
-        f['pix'][~f.mask['pix']]
-
-        Parameters:
-        -----------
-        pixels : If pixels is None, grab all pixels of 'nside_catalog'.
-
-        Returns:
-        --------
-        recarray : pixels and mask value
-        """
-        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:
-            pixels = np.arange(hp.nside2npix(nside_catalog))   
-
-        npix = len(pixels)
-
-        catalog_dir = self['catalog']['dirname']
-        catalog_base = self['catalog']['basename']
-         
-        mask_dir = self['mask']['dirname']
-        mask_base_1 = self['mask']['basename_1']
-        mask_base_2 = self['mask']['basename_2']
-         
-        data = np.ma.empty(npix,dtype=[('pix',int), ('catalog',object), 
-                                       ('mask_1',object), ('mask_2',object)])
-        mask = np.ma.empty(npix,dtype=[('pix',bool), ('catalog',bool), 
-                                       ('mask_1',bool), ('mask_2',bool)])
-        for ii,pix in enumerate(pixels):
-            if pix is None:
-                # DEPRECTATED: ADW 2018-06-17
-                # This is not really being used anymore
-                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)
-            else:
-                catalog = os.path.join(catalog_dir,catalog_base%pix)
-                mask_1 = os.path.join(mask_dir,mask_base_1%pix)
-                mask_2 = os.path.join(mask_dir,mask_base_2%pix)
-            data[ii]['pix'] = pix if pix is not None else -1
-            data[ii]['catalog'] = catalog
-            data[ii]['mask_1']  = mask_1
-            data[ii]['mask_2']  = mask_2
-         
-            mask[ii]['catalog'] = not os.path.exists(catalog)
-            mask[ii]['mask_1']  = not os.path.exists(mask_1)
-            mask[ii]['mask_2']  = not os.path.exists(mask_2)
-
-        for name in ['catalog','mask_1','mask_2']:
-            if np.all(mask[name]): logger.warn("All '%s' files masked"%name)
-
-        # mask 'pix' if all files not present
-        mask['pix'] = mask['catalog'] | mask['mask_1'] | mask['mask_2']
-
-        if np.all(mask['pix']): logger.warn("All pixels masked")
-                
-        #return np.ma.mrecords.MaskedArray(data, mask, fill_value=[-1,None,None,None])
-        #return np.ma.mrecords.MaskedArray(data, mask, fill_value=[-1,'','',''])
-        return np.ma.MaskedArray(data, mask, fill_value=[-1,'','',''])
-
-
     def _createFilenames(self):
         """
         Create a masked records array of all filenames for the given set of
@@ -271,11 +189,17 @@ def _createFilenames(self):
         npix = hp.nside2npix(nside_catalog)
         pixels = np.arange(npix)
 
-        catalog_dir = self['catalog']['dirname']
+        catalog_dir = os.path.expandvars(self['catalog']['dirname'])
+        if not os.path.isdir(catalog_dir):
+            msg = "Directory does not exist: %s"%catalog_dir
+            raise IOError(msg)
         catalog_base = self['catalog']['basename']
         catalog_path = os.path.join(catalog_dir,catalog_base)
 
-        mask_dir    = self['mask']['dirname']
+        mask_dir    = os.path.expandvars(self['mask']['dirname'])
+        if not os.path.isdir(mask_dir):
+            msg = "Directory does not exist: %s"%mask_dir
+            raise IOError(msg)
         mask_base_1 = self['mask']['basename_1']
         mask_base_2 = self['mask']['basename_2']
         mask_path_1 = os.path.join(mask_dir,mask_base_1)
diff --git a/ugali/utils/fileio.py b/ugali/utils/fileio.py
index 9f1caa1..b62bd6b 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.
@@ -105,7 +121,7 @@ def load_header(kwargs):
         if not keys: keys = hdr.keys()
         return {k:hdr[k] for k in keys}
     except Exception as e:
-        logger.error("Failed to load file: %(filename)s"%kwargs)
+        logger.error("Failed to load header from file: %(filename)s"%kwargs)
         raise(e)
 
 def load_headers(filenames,multiproc=False,**kwargs):
@@ -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 a03621d..b4696ce 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?
 
@@ -596,9 +596,9 @@ def read_map(filename, nest=False, hdu=None, h=False, verbose=False):
     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=False):
 
     # 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 7de60e9..5d8cb16 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)
@@ -452,7 +453,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 8ec502f..505530a 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