Fixes for Running at Fermilab

.github/workflows/python-package.yml

@@ -18,7 +18,7 @@ jobs:
     strategy:
       max-parallel: 3
       matrix:
-        python-version: ["3.7", "3.8", "3.9"]
+        python-version: ["3.8", "3.9", "3.10", "3.11"]
     steps:
     - uses: actions/checkout@v4
     - name: Set up Python ${{ matrix.python-version }}

ugali/analysis/farm.py

@@ -128,7 +128,7 @@ def submit(self, pixels, queue=None, debug=False, configfile=None):
         logger.info('=== Submit Likelihood ===')
         for ii,pix in enumerate(pixels):
             msg = '  (%i/%i) pixel=%i nside=%i; (lon, lat) = (%.2f, %.2f)'
-            msg = msg%(ii+1,len(pixels),pix, self.nside_likelihood,lon[ii],lat[ii])
+            msg = msg%(ii,len(pixels),pix, self.nside_likelihood,lon[ii],lat[ii])
             logger.info(msg)
 
             # Create outfile name

ugali/analysis/results.py

@@ -131,6 +131,17 @@ def estimate_position_angle(self,param='position_angle',burn=None,clip=10.0,alph
         return ret
 
     def bayes_factor(self,param,burn=None,clip=10.0,bins=50):
+        """This is a very simplistic calculation of the Bayes Factor from the
+        posterior samples of a single parameter at the specific
+        parameter value of x=0 using the Savage-Dickey method for
+        nested hypothesis. It is currently assumping a flat prior
+        since the prior is not being passed in explicitly.
+
+        More information can be found here:
+        https://www.sciencedirect.com/science/article/abs/pii/S0022249611000666
+
+        ADW 2024-09-06: This should be moved to stats.
+        """
         # CAREFUL: Assumes a flat prior...
         try: 
             data = self.samples.get(param,burn=burn,clip=clip)

ugali/analysis/search.py

@@ -382,11 +382,11 @@ def finalizeObjects(self, objects):
 
         coordsys = self.config['coords']['coordsys']
         if coordsys.lower() == 'gal':
-            print("GAL coordintes")
+            print("GAL coordinates")
             objs['GLON'],objs['GLAT'] = lon,lat
             objs['RA'],objs['DEC'] = gal2cel(lon,lat)
         else:
-            print("CEL coordintes")
+            print("CEL coordinates")
             objs['RA'],objs['DEC'] = lon,lat
             objs['GLON'],objs['GLAT'] = cel2gal(lon,lat)
 

ugali/isochrone/model.py

@@ -219,11 +219,11 @@ def sample(self, mode='data', mass_steps=1000, mass_min=0.1, full_data_range=Fal
             mass_init_min = self.mass_init[self.stage==self.hb_stage].min()
             mass_init_max = self.mass_init[self.stage==self.hb_stage].max()
             cut = (mass_init_array>mass_init_min)&(mass_init_array<mass_init_max)
-            if isinstance(self.hb_spread,collections.Iterable):
+            try: 
                 # Explicit dispersion spacing
                 dispersion_array = self.hb_spread
                 n = len(dispersion_array)
-            else:
+            except TypeError:
                 # Default dispersion spacing
                 dispersion = self.hb_spread
                 spacing = 0.025
@@ -280,8 +280,8 @@ def stellar_mass(self, mass_min=0.1, steps=10000):
 
         d_log_mass = (np.log10(mass_max) - np.log10(mass_min)) / float(steps)
         log_mass = np.linspace(np.log10(mass_min), np.log10(mass_max), steps)
-        mass = 10.**log_mass
-
+        mass = np.clip(10.**log_mass, mass_min, mass_max)
+        
         if mass_min < np.min(self.mass_init):
             mass_act_interpolation = scipy.interpolate.interp1d(np.insert(self.mass_init, 0, mass_min),
                                                                 np.insert(self.mass_act, 0, mass_min))

ugali/observation/roi.py

@@ -152,7 +152,7 @@ def plot(self, value=None, pixel=None):
 
     # ADW: Maybe these should be associated with the PixelRegion objects
     def inPixels(self,lon,lat,pixels):
-        """ Function for testing if coordintes in set of ROI pixels. """
+        """ Function for testing if coordinates in set of ROI pixels. """
         nside = self.config.params['coords']['nside_pixel']
         return ugali.utils.healpix.in_pixels(lon,lat,pixels,nside)
 

ugali/pipeline/run_03.0_likelihood.py

@@ -24,6 +24,8 @@ def run(self):
         farm.submit_all(coords=self.opts.coords,queue=self.opts.queue,debug=self.opts.debug)
 
     if 'merge' in self.opts.run:
+        import numpy as np
+
         logger.info("Running 'merge'...")
         mergefile = self.config.mergefile
         roifile = self.config.roifile
@@ -42,13 +44,13 @@ def run(self):
         for pix in np.unique(superpixel):
             outfile = mergefile%pix
             if exists(outfile) and not self.opts.force:
-                logger.warn("  Found %s; skipping..."%outfile)
+                logger.warning("  Found %s; skipping..."%outfile)
             else:
                 healpix.merge_partial_maps(infiles[superpixel == pix],
                                            outfile,multiproc=8)
 
         if exists(roifile) and not self.opts.force:
-            logger.warn("  Found %s; skipping..."%roifile)
+            logger.warning("  Found %s; skipping..."%roifile)
         else:
             logger.info("  Merging likelihood headers...")
             healpix.merge_likelihood_headers(infiles,roifile)
@@ -64,7 +66,7 @@ def run(self):
         logfile = os.path.join(logdir,'scan_tar.log')
         cmd = 'tar --remove-files -cvzf %s %s'%(tarfile,scanfile)
         if exists(tarfile) and not self.opts.force:
-            logger.warn("  Found %s; skipping..."%tarfile)
+            logger.warning("  Found %s; skipping..."%tarfile)
         else:
             logger.info("  Tarring likelihood files...")
             logger.info(cmd)
@@ -73,17 +75,33 @@ def run(self):
     if 'plot' in self.opts.run:
         # WARNING: Loading the full 3D healpix map is memory intensive.
         logger.info("Running 'plot'...")
+        import numpy as np
         # Should do this in environment variable
         import matplotlib
         matplotlib.use('Agg')
         import pylab as plt
         import ugali.utils.plotting as plotting
-        skymap = ugali.utils.skymap.readSparseHealpixMap(self.config.mergefile,'LOG_LIKELIHOOD')[1]
-        plotting.plotSkymap(skymap)
+        from skymap import DESSkymap
+
+        filenames = self.config.mergefile.split('_%')[0]+'_*.fits'
+        infiles = np.array(sorted(glob.glob(filenames)))
+        hpxmaps = ugali.utils.healpix.read_partial_map(infiles,'LOG_LIKELIHOOD',fullsky=True)[2]
+
         outdir = mkdir(self.config['output']['plotdir'])
-        basename = os.path.basename(self.config.mergefile.replace('.fits','.png'))
-        outfile = os.path.join(outdir,basename)
-        plt.savefig(outfile)
+        basename = os.path.basename(self.config.mergefile.split('_%')[0])
+        for i,hpxmap in enumerate(hpxmaps):
+            #plotting.plotSkymap(hpxmap)
+            smap = DESSkymap()
+            smap.draw_hpxmap(hpxmap, cmap='gray_r', vmax=3.5)
+            smap.draw_inset_colorbar()
+            outfile = os.path.join(outdir,basename+'_%02d.png'%i)
+            print("Writing %s..."%outfile)
+            plt.savefig(outfile)
+        # Make the movie
+        import subprocess
+        cmd = "convert -delay 30 -loop 0 %s/*.png %s.gif"%(outdir,basename)
+        print(cmd)
+        subprocess.check_call(cmd,shell=True)
 
     if 'check' in self.opts.run:
         # Check the completion fraction

ugali/utils/batch.py

@@ -28,7 +28,7 @@
     ('lsf'   ,['express','short','medium','long','xlong','xxl',
                'kipac-ibq','bulletmpi']),
     ('slurm' ,[]),
-    ('condor',['local','vanilla','universe','grid']),
+    ('condor',['local','vanilla']),
 ])
 
 # https://confluence.slac.stanford.edu/x/OaUlCw
@@ -80,9 +80,7 @@ def factory(queue,**kwargs):
     elif name in CLUSTERS['slurm']+QUEUES['slurm']:
         batch = Slurm(**kwargs)
     elif name in CLUSTERS['condor']+QUEUES['condor']:
-        # Need to learn how to use condor first...
         batch = Condor(**kwargs)
-        raise Exception("Condor cluster not implemented")
     else:
         raise TypeError('Unexpected queue name: %s'%name)
 
@@ -105,9 +103,17 @@ def __init__(self, **kwargs):
         self.submit_cmd = "submit %(opts)s %(command)s"
         self.jobs_cmd = "jobs"
 
+    def parse_options(self, **opts):
+        """ Parse command line options. """
+
+        # Default options for the cluster
+        options = odict(self.default_opts)
+        options.update(opts)
+        return ''.join('--%s %s '%(k,v) for k,v in options.items())
+
     def jobs(self):
         out = self.popen(self.jobs_cmd)
-        stdout = out.communicate()[0]
+        stdout = out.communicate()[0].decode()
         return stdout
 
     def njobs(self):
@@ -312,13 +318,30 @@ def parse_options(self, **opts):
 
 class Condor(Batch):
     """ Not implemented yet... """
+
+    _defaults = odict([
+        ('n', '50'),
+    ])
+
+    _mapping = odict([
+        ('jobname','J'),
+        ('logfile','o'),
+        ('njobs','n'),
+    ])
+
     def __init__(self, **kwargs):
         super(Condor,self).__init__(**kwargs)
         logger.warning('Condor cluster is untested')
 
-        self.jobs_cmd = 'condor_q -u %s'%self.username
+        self.jobs_cmd = "cjobs -u %s"%(self.username)
         self.submit_cmd = "csub %(opts)s %(command)s"
 
+    def parse_options(self, **opts):
+        # Default options for the cluster
+        options = odict(self.default_opts)
+        options.update(opts)
+        return ''.join('-%s %s '%(k,v) for k,v in options.items())
+
 if __name__ == "__main__":
     import argparse
     description = __doc__

ugali/utils/fileio.py

@@ -141,7 +141,7 @@ def load_headers(filenames,multiproc=False,**kwargs):
         out = [load_header(kw) for kw in kwargs]
 
     logger.debug('Concatenating headers...')
-    return np.rec.fromrecords([o.values() for o in out],names=out[0].keys())
+    return np.rec.fromrecords([list(o.values()) for o in out], names=list(out[0].keys()))
 
 
 def load_file(kwargs):

ugali/utils/healpix.py

@@ -520,7 +520,7 @@ def merge_likelihood_headers(filenames, outfile, **kwargs):
     data_dict = fileio.load_headers(filenames,ext=ext,keys=keys,mulitproc=8)
     names = data_dict.dtype.names
     data_dict.dtype.names = ['PIXEL' if n=='LKDPIX' else n for n in names]
-        
+
     write_partial_map(outfile, data_dict, nside)
     return data_dict
 

ugali/utils/mlab.py

@@ -17,6 +17,11 @@ def isstring(obj):
     """Python 2/3 compatible string check"""
     return isinstance(obj, six.string_types)
 
+def asscalar(a):
+    """Python 2/3 compatible scalar conversion"""
+    try: return np.asscalar(a)
+    except AttributeError: return a.item()
+
 def rec_append_fields(rec, names, arrs, dtypes=None):
     """
     Return a new record array with field names populated with data

ugali/utils/plotting.py

@@ -24,7 +24,6 @@
 import ugali.utils.config
 import ugali.observation.roi
 import ugali.observation.catalog
-import ugali.utils.skymap
 import ugali.utils.projector
 import ugali.utils.healpix
 import ugali.isochrone
@@ -153,7 +152,7 @@ def sparseHealpixFiles(title, infiles, field='MAGLIM',**kwargs):
     Inputs: field
     """
     #map = ugali.utils.skymap.readSparseHealpixMaps(infiles,field)
-    map = ugali.utils.skymap.read_partial_map(infiles,field)
+    nside,pixels,hpxmap = ugali.utils.healpix.read_partial_map(infiles,field,fullsky=True)
     ax = hp.mollview(map=map, title=title, **kwargs)
     return ax, map
 
@@ -602,7 +601,7 @@ def drawMembership(self, ax=None, radius=None, zidx=0, mc_source_id=1):
 
 
         pix = ang2pix(self.nside, self.lon, self.lat)
-        likelihood_pix = ugali.utils.skymap.superpixel(pix,self.nside,self.config.params['coords']['nside_likelihood'])
+        likelihood_pix = ugali.utils.healpix.superpixel(pix,self.nside,self.config.params['coords']['nside_likelihood'])
         config = self.config
         scan = ugali.analysis.scan.Scan(self.config,likelihood_pix)
         likelihood = scan.likelihood

ugali/utils/projector.py

@@ -8,7 +8,7 @@
 import numpy as np
 
 from ugali.utils.logger import logger
-from ugali.utils.mlab import isstring
+from ugali.utils.mlab import isstring, asscalar
 
 ############################################################
 
@@ -117,7 +117,7 @@ def sphereToImage(self, lon, lat):
         lon_rotated, lat_rotated = self.rotator.rotate(lon.flat, lat.flat)
         x, y = self.sphere_to_image_func(lon_rotated, lat_rotated)
 
-        if scalar: return np.asscalar(x), np.asscalar(y)
+        if scalar: return asscalar(x), asscalar(y)
         else:      return x.reshape(lon.shape), y.reshape(lat.shape)
 
     sphere2image = sphereToImage
@@ -129,7 +129,7 @@ def imageToSphere(self, x, y):
         lon_rotated, lat_rotated = self.image_to_sphere_func(x.flat, y.flat)
         lon, lat = self.rotator.rotate(lon_rotated, lat_rotated, invert = True)
 
-        if scalar: return np.asscalar(lon), np.asscalar(lat)
+        if scalar: return asscalar(lon), asscalar(lat)
         else:      return lon.reshape(x.shape), lat.reshape(y.shape)
 
     image2sphere = imageToSphere

ugali/utils/skymap.py

@@ -6,7 +6,8 @@
 import healpy as hp
 
 import ugali.utils.projector
-from ugali.utils.healpix import superpixel,subpixel,ang2pix,pix2ang,query_disc
+from ugali.utils.healpix import superpixel,subpixel
+from ugali.utils.healpix import ang2pix,pix2ang,query_disc
 from ugali.utils.healpix import read_partial_map
 from ugali.utils.logger import logger
 from ugali.utils.config import Config