Merge pull request DarkEnergySurvey#90 from kadrlica/y3-mw-dev

Updates from MW Satellite Census 2020

tests/test_results.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+"""
+Generic python script.
+"""
+__author__ = "Alex Drlica-Wagner"
+
+import numpy as np
+import ugali.analysis.results
+
+def test_surface_brightness():
+
+    # Draco values from McConnachie et al. 2012
+    abs_mag = -8.8
+    app_mag = 10.6
+    distance = 76 #kpc
+    distance_modulus = 19.40
+    a_half_arcmin = 10.0 # arcmin
+    a_physical_kpc = 0.221 # kpc
+    ellip = 0.31
+    mu0 = 26.1 # mag/arcsec^2
+
+    # Convert to azimuthally average
+    r_half_arcmin = a_half_arcmin * np.sqrt(1-ellip)
+    r_physical_kpc = a_physical_kpc * np.sqrt(1-ellip)
+
+    mu1 = ugali.analysis.results.surfaceBrightness(abs_mag, r_physical_kpc, distance)
+    mu2 = ugali.analysis.results.surfaceBrightness2(app_mag, r_half_arcmin)
+
+    np.testing.assert_allclose(mu1,mu0,atol=2e-2)
+    np.testing.assert_allclose(mu2,mu0,atol=2e-2)

ugali/analysis/farm.py

@@ -54,52 +54,6 @@ def command(self, outfile, configfile, pix):
                       verbose='-v' if self.verbose else '')
         cmd = '%(script)s %(config)s %(outfile)s --hpx %(nside)i %(pix)i %(verbose)s'%params
         return cmd
-
-    # ADW: Should probably be in a utility
-    def footprint(self, nside=None):
-        """
-        UNTESTED.
-        Should return a boolean array representing the pixels in the footprint.
-        """
-        if nside is None:
-            nside = self.nside_pixel
-        elif nside < self.nside_catalog: 
-            raise Exception('Requested nside=%i is less than catalog_nside'%nside)
-        elif nside > self.nside_pixel:
-            raise Exception('Requested nside=%i is greater than pixel_nside'%nside)
-        pix = np.arange(hp.nside2npix(nside), dtype=int)
-        map = self.inFootprint(pix,nside)
-        return map 
-
-    # ADW: Should probably be in a utility
-    def inFootprint(self, pixels, nside=None):
-        """
-        Open each valid filename for the set of pixels and determine the set 
-        of subpixels with valid data.
-        """
-        if np.isscalar(pixels): pixels = np.array([pixels])
-        if nside is None: nside = self.nside_likelihood
-
-        inside = np.zeros( len(pixels), dtype='bool')
-        if not self.nside_catalog:
-            catalog_pix = [0]
-        else:
-            catalog_pix = superpixel(pixels,nside,self.nside_catalog)
-            catalog_pix = np.intersect1d(catalog_pix,self.catalog_pixels)
-
-        for filenames in self.filenames[catalog_pix]:
-            # ADW: Need to replace with healpix functions...
-            #logger.debug("Loading %s"%filenames['mask_1'])
-            #subpix_1,val_1 = ugali.utils.skymap.readSparseHealpixMap(filenames['mask_1'],'MAGLIM',construct_map=False)
-            _n,subpix_1,val_1 = read_partial_map(filenames['mask_1'],'MAGLIM',fullsky=False)
-            #logger.debug("Loading %s"%filenames['mask_2'])
-            #subpix_2,val_2 = ugali.utils.skymap.readSparseHealpixMap(filenames['mask_2'],'MAGLIM',construct_map=False)
-            _n,subpix_2,val_2 = read_partial_map(filenames['mask_2'],'MAGLIM',fullsky=False)
-            subpix = np.intersect1d(subpix_1,subpix_2)
-            superpix = np.unique(superpixel(subpix,self.nside_pixel,nside))
-            inside |= np.in1d(pixels, superpix)
-
-        return inside
 
     def submit_all(self, coords=None, queue=None, debug=False):
         """
@@ -114,18 +68,9 @@ def submit_all(self, coords=None, queue=None, debug=False):
         if coords is None:
             pixels = np.arange(hp.nside2npix(self.nside_likelihood))
         else:
-            coords = np.asarray(coords)
-            if coords.ndim == 1:
-                coords = np.array([coords])
-            if coords.shape[1] == 2:
-                lon,lat = coords.T
-                radius  = np.zeros(len(lon))
-            elif coords.shape[1] == 3:
-                lon,lat,radius = coords.T
-            else:
-                raise Exception("Unrecognized coords shape:"+str(coords.shape))
+            lon,lat,radius = coords['lon'],coords['lat'],coords['radius']
 
-            #ADW: targets is still in glon,glat
+            #ADW: coords are always parsed in GAL, so convert to CEL if necessary
             if self.config['coords']['coordsys'].lower() == 'cel':
                 lon,lat = gal2cel(lon,lat)
 
@@ -158,11 +103,11 @@ def submit(self, pixels, queue=None, debug=False, configfile=None):
         """
         # For backwards compatibility
         batch = self.config['scan'].get('batch',self.config['batch'])
-        queue = batch['cluster'] if queue is None else queue
+        queue = batch.get('default','medium') if queue is None else queue
 
         # Need to develop some way to take command line arguments...
-        self.batch = ugali.utils.batch.batchFactory(queue,**batch['opts'])
-        self.batch.max_jobs = batch.get('max_jobs',200)
+        self.batch = ugali.utils.batch.batchFactory(queue,**batch.get(queue,{}))
+        self.batch.max_jobs = self.config['scan'].get('max_jobs',200)
 
         if np.isscalar(pixels): pixels = np.array([pixels])
 
@@ -178,7 +123,7 @@ def submit(self, pixels, queue=None, debug=False, configfile=None):
 
         lon,lat = pix2ang(self.nside_likelihood,pixels)
         commands = []
-        chunk = batch['chunk']
+        chunk = self.config['scan'].get('chunk',25)
         istart = 0
         logger.info('=== Submit Likelihood ===')
         for ii,pix in enumerate(pixels):
@@ -189,7 +134,7 @@ def submit(self, pixels, queue=None, debug=False, configfile=None):
             # Create outfile name
             outfile = self.config.likefile%(pix,self.config['coords']['coordsys'].lower())
             outbase = os.path.basename(outfile)
-            jobname = batch['jobname']
+            jobname = batch.get('jobname','ugali')
 
             # Submission command
             sub = not os.path.exists(outfile)
@@ -226,21 +171,48 @@ def submit(self, pixels, queue=None, debug=False, configfile=None):
                 time.sleep(0.5)
 
     def write_script(self, filename, commands):
+        """ Write a batch submission script.
+
+        Parameters
+        ----------
+        filename : filename of batch script
+        commands : list of commands to execute
+
+        Returns
+        -------
+        None
+        """
         info = 'echo "{0:=^60}";\n'
         hline = info.format("")
         newline = 'echo;\n'
+        shebang = "#!/usr/bin/env bash"
+        # Limit the memory based on SLAC 4 GB per node (defined in KB)
+        # Careful, shell arithmetic is weird.
+        memory_limit = """
+if [ -n "$LSB_CG_MEMLIMIT" ] & [ -n "$LSB_HOSTS" ]; then
+    mlimit=$(( $(wc -w <<< $LSB_HOSTS) * $LSB_CG_MEMLIMIT/1024 * 9/10 ))
+    ulimit -v ${mlimit}; ulimit -H -v ${mlimit};
+fi
+"""
+        memory_usage=r"""free -m | awk 'NR==2{printf "Memory Usage: %.2f/%.2fGB (%.2f%%)\n",$3/1024,$2/1024,$3*100/$2}';"""
+        memory_usage=r"""ps -U $USER --no-headers -o rss | awk '{sum+=$1} END {print "Memory Usage: " int(sum/1024**2) "GB"}'"""
 
         istart, iend = commands[0][0], commands[-1][0]
         script = open(filename,'w')
+        script.write(shebang)
+        #script.write(memory_limit)
         script.write(hline)
-        script.write(info.format('Submit Pixels %i to %i'%(istart,iend)))
+        script.write(info.format('Submit Jobs %i to %i'%(istart,iend)))
         script.write(hline)
         script.write(newline)
         script.write('status=0;\n')
         for i,cmd,lon,lat,sub in commands: 
-            script.write(info.format('Pixel %i: (%.2f, %.2f)'%(i,lon,lat)))
-            if sub: script.write('%s; [ $? -ne 0 ] && status=1;\n'%cmd)
-            else:   script.write('echo "%s";\n'%self.skip)
+            script.write(info.format('Job %i: (%.2f, %.2f)'%(i,lon,lat)))
+            if sub: 
+                script.write(memory_usage+'\n')
+                script.write('%s; [ $? -ne 0 ] && status=1;\n'%cmd)
+            else:   
+                script.write('echo "%s";\n'%self.skip)
             script.write(hline)
             script.write(newline)
         script.write('exit $status;\n')

ugali/analysis/kernel.py

@@ -40,29 +40,42 @@ def __init__(self, proj='ait', **kwargs):
         super(Kernel,self).__init__(**kwargs)
 
     def __call__(self, lon, lat):
+        """
+        Return the value of the pdf at a give location.
+
+        Parameters
+        ----------
+        lon : longitude (deg)
+        lat : latitude (deg)
+        
+        Returns
+        -------
+        pdf : normalized, truncated pdf
+        """
         return self.pdf(lon, lat)
 
     @abstractmethod
     def _kernel(self, r):
-        # unnormalized, untruncated kernel
+        """Unnormalized, untruncated kernel"""
         pass
 
     def _pdf(self, radius):
-        # unnormalized, truncated kernel
+        """Unnormalized, truncated kernel"""
         return np.where(radius<=self.edge, self._kernel(radius), 0.)
 
     @abstractmethod
     def pdf(self, lon, lat):
-        # normalized, truncated pdf
+        """Normalized, truncated pdf"""
         pass
 
     @property
     def norm(self):
-        # Numerically integrate the pdf
+        """Normalization from the integated pdf"""
         return 1./self.integrate()
 
     @property
     def projector(self):
+        """Projector used to transform to local sky coordinates."""       
         if self.proj is None or self.proj.lower()=='none':
             return None
         else:
@@ -133,12 +146,14 @@ class EllipticalKernel(Kernel):
     _params = odict(
         list(Kernel._params.items()) + 
         [
+            # This is the semi-major axis
             ('extension',     Parameter(0.1, [0.0001,0.5]) ),   
-            ('ellipticity',   Parameter(0.0, [0.0, 0.99]) ),    # Default 0 for RadialKernel
-            ('position_angle',Parameter(0.0, [0.0, 180.0]) ),  # Default 0 for RadialKernel
-            # This is the PA *WEST* of North.
+            # This is e = 1 - b/a (0 for RadialKernel)
+            ('ellipticity',   Parameter(0.0, [0.0, 0.99]) ),    
+            # This is the PA *WEST* of North (0 for RadialKernel)
             # to get the conventional PA EAST of North take 90-PA
             # Would it be better to have bounds [-90,90]?
+            ('position_angle',Parameter(0.0, [0.0, 180.0]) ),  # 
         ])
     _mapping = odict([
         ('e','ellipticity'),
@@ -183,19 +198,36 @@ def sample_radius(self, n):
         """
         Sample the radial distribution (deg) from the 2D stellar density.
         Output is elliptical radius in true projected coordinates.
+
+        Parameters
+        ----------
+        n : number of stars to sample
+
+        Returns
+        -------
+        radius : distance from centroid (deg)
         """
+        size = int(n)
         edge = self.edge if self.edge<20*self.extension else 20*self.extension
-        radius = np.linspace(0, edge, 1.e5)
+        radius = np.linspace(0, edge, int(10**5))
         pdf = self._pdf(radius) * np.sin(np.radians(radius))
         cdf = np.cumsum(pdf)
         cdf /= cdf[-1]
         fn = scipy.interpolate.interp1d(cdf, list(range(0, len(cdf))))
-        index = np.floor(fn(np.random.uniform(size=n))).astype(int)
+        index = np.floor(fn(np.random.uniform(size=size))).astype(int)
         return radius[index]
 
     def sample_lonlat(self, n):
         """
         Sample 2D distribution of points in lon, lat
+
+        Parameters
+        ----------
+        n : number of stars to sample
+
+        Returns
+        -------
+        lon,lat : longitude and latitude (deg)
         """
         # From http://en.wikipedia.org/wiki/Ellipse#General_parametric_form
         # However, Martin et al. (2009) use PA theta "from North to East"
@@ -374,7 +406,7 @@ def _kernel(self, radius):
         return ((1./np.sqrt(1.+(radius/self.r_c)**2))-(1./np.sqrt(1.+(self.r_t/self.r_c)**2)))**2
 
     def _cache(self, name=None):
-        if name in ['extension','ellipticity','truncate']:
+        if name in [None,'extension','ellipticity','truncate']:
             self._norm = 1./self.integrate()
         else:
             return

ugali/analysis/model.py

@@ -25,9 +25,9 @@ def asscalar(a):
     # Do we want to check that the value is numeric?
     #if   isinstance(value, (int, long, float)): return value
     try:
-        return np.asscalar(a)
+        return a.item()
     except AttributeError as e:
-        return np.asscalar(np.asarray(a))
+        return np.asarray(a).item()
 
 
 class Model(object):