Merge lp:~maddm/maddm/no_numpy_scipy into lp:maddm

[Olivier Mattelaer]

Avoid direct crash if numpy/scipy are not present on the system.
-> forbids indirect detection but allows relic density computation and direct detection.

This allows to have maddm v2 options still working in version 3 even if numpy/scipy are not available.

[lp:~maddm/maddm/no_numpy_scipy]

18. By olivier-mattelaer on 2018-05-06

some issue created by the change

19. By olivier-mattelaer on 2018-05-06

additional fix

20. By Federico Ambrogi <email address hidden> on 2018-05-15

Formatted string as as python 2.6

21. By olivier-mattelaer on 2018-05-16

change according to federico review

22. By olivier-mattelaer on 2018-05-16

remove some printout

Unmerged revisions

=== modified file '__init__.py'
--- __init__.py	2018-03-30 22:19:49 +0000
+++ __init__.py	2018-05-06 15:24:43 +0000
@@ -16,6 +16,6 @@
 new_reweight = {'indirect': MGoutput.Indirect_Reweight}
 
 ## The test/code have been validated up to this version
-latest_validated_version = (2,6,2)
+latest_validated_version = (2,6,3)
 minimal_mg5amcnlo_version = (2,6,2)
 maximal_mg5amcnlo_version = (1000,1000,1000)

=== modified file 'maddm_interface.py'
--- maddm_interface.py	2018-04-29 20:55:19 +0000
+++ maddm_interface.py	2018-05-06 15:24:43 +0000
@@ -115,6 +115,11 @@
     def post_install_PPPC4DMID(self):
         if os.path.exists(pjoin(MG5DIR, 'PPPC4DMID')):
             self.options['pppc4dmid_path'] = pjoin(MG5DIR, 'PPPC4DMID')
+        
+        if not maddm_run_interface.HAS_SCIPY:
+            logger.critical("PPC4DMID module requires scipy to be working. Please install those python module. (they are not present)")
+            logger.info("you can try to use \"pip install scipy\"")
+            
         return
     
     def set_configuration(self, config_path=None, final=True, **opts):

=== modified file 'maddm_run_interface.py'
--- maddm_run_interface.py	2018-04-29 20:53:38 +0000
+++ maddm_run_interface.py	2018-05-06 15:24:43 +0000
@@ -10,12 +10,7 @@
 import stat
 import shutil
 
-from scipy.interpolate import interp1d
-from scipy.integrate import quad
-from scipy.optimize import minimize_scalar
-from scipy.optimize import brute
-from scipy.optimize import fmin
-from scipy.special import gammainc
+
 
 import MGoutput
 from madgraph import MadGraph5Error
@@ -34,18 +29,36 @@
         
 #import darkmatter as darkmatter
 
-import numpy as np
+logger = logging.getLogger('madgraph.plugin.maddm')
 
 try:
     import pymultinest
-except:
-    print('WARNING: Multinest module not found! All multinest parameter scanning features will be disabled.')
-
-
-#import types
+except ImportError:
+    pass
+
+try:
+    from scipy.interpolate import interp1d
+    from scipy.integrate import quad
+    from scipy.special import gammainc
+except ImportError, error:
+    print error
+    logger.warning('scipy module not found! Some Indirect detection features will be disabled.')
+    HAS_SCIPY = False 
+else:
+    HAS_SCIPY = True
+
+try:
+    import numpy as np    
+except ImportError:
+    logger.warning('numpy module not found! Indirect detection features will be disabled.')
+    HAS_NUMPY = False
+else:
+    HAS_NUMPY = True
+        
+class ModuleMissing(Exception): pass
 
 pjoin = os.path.join
-logger = logging.getLogger('madgraph.plugin.maddm')
+
 logger_tuto = logging.getLogger('tutorial_plugin')
 #logger.setLevel(10) #level 20 = INFO
 
@@ -111,14 +124,18 @@
             self._sigma_ID[item] = -1.0
             self._sigma_ID_width[item] = -1.0
 
-        self.load_constraints()
+        if HAS_NUMPY:
+            self.load_constraints()
+
 
         logger.info('Loaded experimental constraints. To change, use the set command')
         logger.info('Omega h^2 = %.4e +- %.4e' %(self._oh2_planck, self._oh2_planck_width))
-        logger.info('Spin Independent cross section: %s' % self._dd_si_limit_file)
-        logger.info('Spin Dependent cross section (p): %s' % self._dd_sd_proton_limit_file)
-        logger.info('Spin Dependent cross section (n): %s' % self._dd_sd_neutron_limit_file)
-
+        if HAS_NUMPY:
+            logger.info('Spin Independent cross section: %s' % self._dd_si_limit_file)
+            logger.info('Spin Dependent cross section (p): %s' % self._dd_sd_proton_limit_file)
+            logger.info('Spin Dependent cross section (n): %s' % self._dd_sd_neutron_limit_file)
+        else:
+            logger.info('Spin (in)dependent not available due to the missing python module: numpy')
 #        for chan in self._allowed_final_states:
 #            logger.info('Indirect Detection cross section for final state %s at velocity %.2e: %s'\
 #                        % (chan, self._id_limit_vel[chan] ,self._id_limit_file[chan]))
@@ -146,6 +163,10 @@
 
     #Returns a value in cm^2
     def SI_max(self, mdm):
+        if not HAS_NUMPY:
+            logger.warning("missing numpy module for SI limit")
+            return __infty__
+        
         if (mdm < np.min(self._dd_si_limit_mDM) or mdm > np.max(self._dd_si_limit_mDM)):
             logger.warning('Dark matter mass value '+str(mdm)+' is outside the range of SI limit')
             return __infty__
@@ -154,6 +175,10 @@
 
     #Returns a value in cm^2
     def SD_max(self,mdm, nucleon):
+        if not HAS_NUMPY:
+            logger.warning("missing numpy module for SD limit")
+            return __infty__
+            
         if nucleon not in ['n','p']:
             logger.error('nucleon can only be p or n')
             return __infty__
@@ -173,6 +198,9 @@
 
     #Returns a value in cm^3/s
     def ID_max(self,mdm, channel):
+        if not HAS_NUMPY:
+            logger.warning("missing numpy module for ID limit")
+            return __infty__
         if (mdm < np.min(self._id_limit_mdm[channel]) or mdm > np.max(self._id_limit_mdm[channel])):
             logger.warning('Dark matter mass value %.2e for channel %s is outside the range of ID limit' % (mdm, channel))
             return __infty__
@@ -267,6 +295,10 @@
     # this function extracts the values of the spectra interpolated linearly between two values mdm_1 and mdm_2                                              
     # mdm is the DM candidate mass, spectrum is gammas, positron etc, channel is the SM annihilation e.g. bbar, hh etc.                        
     def interpolate_spectra(self, sp_dic, mdm = '', spectrum = '' , channel = '' , earth = False , prof = 'Ein' , prop = 'MED' , halo_func = 'MF1'):
+        
+        if not HAS_SCIPY:
+            raise ModuleMissing('scipy module is required for this functionality.')
+        
         M   = sp_dic['Masses']
         dm_min =  max([m for m in M if m <= mdm])  # extracting lower mass limit to interpolate from                                                                  
         dm_max =  min([m for m in M if m >= mdm])  # extracting upper mass limit to interpolate from                                                                    
@@ -306,7 +338,8 @@
         self.dSph_ll_files_path =  pjoin(MDMDIR, 'Fermi_Data', 'likelihoods')
         self.dwarves_list = ['coma_berenices', 'draco', 'segue_1', 'ursa_major_II', 'ursa_minor', 'reticulum_II' ] # dSphs with the 6 highest Jfactors
         self.dwarveslist_all = self.extract_dwarveslist() 
-        self.dw_in = self.dw_dic()
+        if HAS_NUMPY:
+            self.dw_in = self.dw_dic()
         self.ll_tot = ''
 
     # This function reads the list of dwarves from the Jfactor.dat file
@@ -363,6 +396,10 @@
     def eflux(self,spectrum, emin=1e2, emax=1e5, quiet=False):
         """ Integrate a generic spectrum, multiplied by E, to get the energy flux.                                                                                             
         """
+        
+        if not HAS_SCIPY:
+            raise ModuleMissing('scipy module is required for this functionality.')
+        
         espectrum = lambda e: spectrum(e)*e
         tol = min(espectrum(emin),espectrum(emax))*1e-10
         try:
@@ -373,6 +410,9 @@
 
     def marg_like_dw(self,dw_in_i,pred,marginalize):
       
+        if not HAS_SCIPY:
+            raise ModuleMissing('scipy module is required for this functionality.')
+        
         j0, nBin = self.j0 , self.nBin
 
         j,jerr,like_inter = dw_in_i['Jfac'], dw_in_i['Jfac_err'], dw_in_i['likelihood']
@@ -403,6 +443,9 @@
 
     def res_tot_dw(self,pred,marginalize):
 
+        if not HAS_SCIPY:
+            raise ModuleMissing('scipy module is required for this functionality.')
+
         dw_in = self.dw_in
         ll_tot = 0.0
         ll_null = 0.0
@@ -437,6 +480,13 @@
     def Fermi_sigmav_lim(self, mDM, x = '' , dndlogx = '' , marginalize = True, sigmav_th = False , maj_dirac='', \
                                sigmavmin=1e-35, sigmavmax=1e-15, step_size_scaling=1.0, cl_val = 0.95):
 
+        if not HAS_NUMPY:
+            logger.warning("Fermi limit ignored due to missing numpy module")
+            return -1
+        if not HAS_SCIPY:
+            logger.warning("Fermi limit ignored due to missing scipy module")
+            return -1        
+        
         np.seterr(divide='ignore', invalid='ignore')   # Keep numpy from complaining about dN/dE = 0...                                                                
         j0 , nBin = self.j0 , self.nBin
         
@@ -1226,8 +1276,8 @@
             sp_name = sp + '_spectrum_pythia8.dat'
             out_dir = pjoin(self.dir_path,'Indirect', 'Events', run_name, sp_name )
             if sp == 'gammas': # x values are the same for all the spectra
-                  x = np.loadtxt(out_dir , unpack = True )[0]
-                  self.Spectra.spectra['x'] = [ np.power(10,num) for num in x]     # from log[10,x] to x
+                x = np.loadtxt(out_dir , unpack = True )[0]
+                self.Spectra.spectra['x'] = [ np.power(10,num) for num in x]     # from log[10,x] to x
             self.Spectra.spectra[sp] = np.loadtxt(out_dir , unpack = True )[1].tolist()                                  
        
 
@@ -1239,6 +1289,10 @@
             logger.error('PPPC4DMID not installed. Please install by typing "install PPPC4DMID".')
             return
 
+        if not HAS_SCIPY:
+            logger.error('using PPPC4DMID requires scipy module. Please install it (for example with "pip install scipy")')
+            return
+
         if 'PPPC4DMID' in self.maddm_card['indirect_flux_source_method'] or 'inclusive' in self.maddm_card['sigmav_method']:
             if self.Spectra.check_mass(mdm):
                if '_ew' in self.maddm_card['indirect_flux_source_method']:
@@ -2356,7 +2410,9 @@
     def set_default_indirect(self):
         """set the default value for relic="""
         
-        if self.availmode['has_indirect_detection']:
+        if not HAS_NUMPY:
+            self.switch['indirect'] = 'Not Avail. (numpy missing)'
+        elif self.availmode['has_indirect_detection']:
             self.switch['indirect'] = 'sigmav'     
         else:
             self.switch['indirect'] = 'Not Avail.'
@@ -2368,7 +2424,9 @@
         if hasattr(self, 'allowed_indirect'):
             return getattr(self, 'allowed_indirect')
 
-        if self.availmode['has_indirect_detection']:
+        if not HAS_NUMPY:
+            self.allowed_indirect =  ['OFF']
+        elif self.availmode['has_indirect_detection']:
             self.allowed_indirect =  ['OFF', 'sigmav', 'flux_source', 'flux_earth']
         else:
             return []