Merge pull request #11 from fermiPy/dmcat-dev

commits for glory duck

dmpipe/dm_prepare.py

@@ -9,6 +9,7 @@
 from __future__ import absolute_import, division, print_function
 
 import os
+import sys
 import copy
 
 from shutil import copyfile
@@ -167,6 +168,13 @@ def _write_profile_yaml(cls, target, profile_path, targ_ver, spatial):
             source_model.update(dict(SpatialModel='DiffuseSource',
                                      SpatialType='SpatialMap',
                                      Spatial_Filename=target.j_map_file))
+        elif spatial in ['hpxmap']:
+            if target.j_map_file is None:
+                j_map_file = profile_path.replace('.yaml', '.fits')
+                target.write_jmap_hpx(j_map_file)
+            source_model.update(dict(SpatialModel='DiffuseSource',
+                                     SpatialType='SpatialMap',
+                                     Spatial_Filename=target.j_map_file))
         elif spatial in ['radial']:
             target.j_rad_file = profile_path.replace('.yaml', '.dat')
             target.write_j_rad_file()
@@ -180,6 +188,13 @@ def _write_profile_yaml(cls, target, profile_path, targ_ver, spatial):
             source_model.update(dict(SpatialModel='DiffuseSource',
                                      SpatialType='SpatialMap',
                                      Spatial_Filename=target.d_map_file))
+        elif spatial in ['hpxdmap']:
+            if target.d_map_file is None:
+                d_map_file = profile_path.replace('.yaml', '.fits')
+                target.write_dmap_hpx(d_map_file)
+            source_model.update(dict(SpatialModel='DiffuseSource',
+                                     SpatialType='SpatialMap',
+                                     Spatial_Filename=target.d_map_file))
         elif spatial in ['dradial']:
             target.d_rad_file = profile_path.replace('.yaml', '.dat')
             target.write_d_rad_file()

dmpipe/dm_spectral.py

@@ -60,12 +60,12 @@ class ConvertCastro(Link):
     @staticmethod
     def is_decay_sed(sedfile):
         tokens = os.path.splitext(os.path.basename(sedfile))[0].split('_')
-        return tokens[2] in ['point', 'dmap', 'dradial']
+        return tokens[2] in ['point', 'dmap', 'dradial', 'hpxdmap']
 
     @staticmethod
     def is_ann_sed(sedfile):
         tokens = os.path.splitext(os.path.basename(sedfile))[0].split('_')
-        return tokens[2] in ['point', 'map', 'radial']
+        return tokens[2] in ['point', 'map', 'radial', 'hpxmap']
 
     @staticmethod
     def select_channels(channels, sedfile):

dmpipe/dm_spectral_utils.py

@@ -696,6 +696,114 @@ def x_edges(self):
         return np.power(10, log_masses)
 
 
+    def write_glory_duck_v1(self, basepath, **kwargs):
+        """Write this object into the format that the 'Glory Duck' collaboration wants
+
+        Paramters
+        ---------
+
+        basepath : str
+            Path to output files
+        """
+        kwcopy = kwargs.copy()
+        sigmav_min = kwcopy.get('logsigmav_min', -30)
+        sigmav_max = kwcopy.get('logsigmav_max', -24)
+        #sigmav_nstep = kwcopy.get('sigmav_nstep', 20000)
+        sigmav_nstep = kwcopy.get('sigmav_nstep', 200)
+
+        sigmav_steps = np.logspace(sigmav_min, sigmav_max, sigmav_nstep)[::-1]
+
+        for imass, mass in enumerate(self.masses):
+            filepath = "%s_%0.1f_GeV.txt" % (basepath, mass)
+
+            ll_interp = self[imass].interp
+            ll_vals = ll_interp(sigmav_steps)
+
+            outfile = open(filepath, 'w!')
+
+            sys.stdout.write("Writing %i values to %s: " % (sigmav_nstep, filepath))
+            sys.stdout.flush()
+            for i, (sigmav, ll) in enumerate(zip(sigmav_steps, ll_vals)):
+                if i % 100 == 0:
+                    sys.stdout.write('.')
+                    sys.stdout.flush()
+                outfile.write("%6e %6e\n" % (sigmav, ll))
+            outfile.close()
+            sys.stdout.write('!\n')
+
+
+    def write_glory_duck_v2(self, filepath, **kwargs):
+        """Write this object into the format that the 'Glory Duck' collaboration wants
+
+        Paramters
+        ---------
+
+        basepath : str
+            Path to output files
+        """
+        kwcopy = kwargs.copy()
+        sigmav_min = kwcopy.get('logsigmav_min', -30)
+        sigmav_max = kwcopy.get('logsigmav_max', -24)
+        #sigmav_nstep = kwcopy.get('sigmav_nstep', 20000)
+        sigmav_nstep = kwcopy.get('sigmav_nstep', 200)
+        mass_vals = kwcopy.get('mass_vals', None)
+
+        sigmav_steps = np.logspace(sigmav_min, sigmav_max, sigmav_nstep)[::-1]
+
+        ll_vals_list = []
+        if mass_vals is not None:
+            mass_list = mass_vals
+        else:
+            mass_list = self.masses
+
+        # Do linear interpolation in log-space
+        log_masses = np.log10(self.masses)
+        log_mass_diffs = log_masses[1:] - log_masses[0:-1]
+        log_mass_vals = np.log10(mass_list)        
+
+        for log_mass in log_mass_vals:
+            mass_index = np.searchsorted(log_masses, log_mass)
+            if mass_index < 0: 
+                ll_vals_list.append(np.zeros((sigmav_nstep)))
+                continue
+            elif mass_index >= len(log_masses):
+                ll_vals_list.append(np.zeros((sigmav_nstep)))
+                continue
+            elif mass_index == len(log_masses) - 1:
+                ll_interp = self[mass_index].interp
+                ll_vals = ll_interp(sigmav_steps)
+                ll_vals_list.append(ll_vals)
+                continue
+            frac_lo = (log_masses[mass_index+1] - log_mass) / log_mass_diffs[mass_index]
+            frac_hi = 1 - frac_lo
+            ll_interp_lo = self[mass_index].interp            
+            ll_vals_lo = ll_interp_lo(sigmav_steps)
+            ll_interp_hi = self[mass_index+1].interp            
+            ll_vals_hi = ll_interp_hi(sigmav_steps)            
+            ll_vals = frac_lo*ll_vals_lo + frac_hi*ll_vals_hi
+            ll_vals_list.append(ll_vals)
+
+        ll_vals_stack = np.vstack(ll_vals_list).T
+
+        outfile = open(filepath, 'w!')      
+        sys.stdout.write("Writing %i values to %s: " % (sigmav_nstep, filepath))
+
+        outfile.write("%.2f" % np.log10(self._astro_value))
+        for mass in mass_list:
+            outfile.write(" %.2f" % (mass))
+        outfile.write("\n")
+        sys.stdout.flush()
+        for i, (sigmav, ll_vals) in enumerate(zip(sigmav_steps, ll_vals_stack)):
+            if i % 1000 == 0:
+                sys.stdout.write('.')
+                sys.stdout.flush()
+            outfile.write("%6e " % (sigmav))
+            for ll_val in ll_vals:
+                outfile.write("%6e " % (ll_val))
+            outfile.write("\n")
+        outfile.close()
+        sys.stdout.write('!\n')
+
 
 class DMSpecTable(object):
     """ Version of the DM spectral tables in tabular form
@@ -857,6 +965,8 @@ def ebin_refs(self):
         """
         return self._e_table["E_REF"].data
 
+
+
     def write_fits(self, filepath, clobber=False):
         """ Write this object to a FITS file