Merge pull request #31 from mirochaj/dm_extension_branch

Functionality added includes interpolation between existing WDM HMF tables, option to use MAR inferred from CDM even when WDM being used.

README.md

@@ -143,3 +143,4 @@ Additional contributions / corrections / suggestions from:
 - Matteo Leo
 - Emma Klemets
 - Felix Bilodeau-Chagnon
+- Joshua Hibbard

ares/physics/HaloMassFunction.py

@@ -55,16 +55,17 @@
     from hmf import MassFunction
     have_hmf = True
     hmf_vstr = hmf.__version__
-    hmf_vers = float(hmf_vstr[0:hmf_vstr.rfind('.')])
+    hmf_vers = float(hmf_vstr[0:hmf_vstr.index('.')+2])
 except ImportError:
     have_hmf = False
     hmf_vers = 0
 
-if 0 < hmf_vers < 3.1:
-    try:
-        from hmf.wdm import MassFunctionWDM
-    except ImportError:
-        pass
+if have_hmf:
+	if 0 <= hmf_vers <= 3.4:
+		try:
+		    MassFunctionWDM = hmf.wdm.MassFunctionWDM
+		except ImportError:
+		    pass
 
 # Old versions of HMF
 try:
@@ -83,6 +84,7 @@
 
 sqrt2 = np.sqrt(2.)
 
+tiny_dndm = 1e-30
 tiny_fcoll = 1e-18
 tiny_dfcolldz = 1e-18
 
@@ -155,6 +157,7 @@ def __init__(self, **kwargs):
             _path = '{0!s}/input/hmf'.format(ARES)
 
         # Look for tables in input directory
+
         if ARES is not None and self.pf['hmf_load'] and (self.tab_name is None):
             prefix = self.tab_prefix_hmf(True)
             fn = '{0!s}/{1!s}'.format(_path, prefix)
@@ -171,6 +174,7 @@ def __init__(self, **kwargs):
                 prefix = self.tab_prefix_hmf()
                 candidates =\
                     glob.glob('{0!s}/input/hmf/{1!s}*'.format(ARES, prefix))
+                print(candidates)
 
                 if len(candidates) == 1:
                     self.tab_name = candidates[0]
@@ -215,6 +219,7 @@ def __init__(self, **kwargs):
 
                         self.tab_name = candidate
 
+
         # Override switch: compute Press-Schechter function analytically
         if self.hmf_func == 'PS' and self.hmf_analytic:
             self.tab_name = None
@@ -282,9 +287,39 @@ def __getattr__(self, name):
         return self.__dict__[name]
 
     def _load_hmf_wdm(self): # pragma: no cover
+
         m_X = self.pf['hmf_wdm_mass']
 
-        _fn = self.tab_prefix_hmf(True)
+        if self.pf['hmf_wdm_interp']:
+            wdm_file_hmfs = []
+            import glob
+            for wdm_file in glob.glob('{!s}/input/hmf/*'.format(ARES)):
+                if self.pf['hmf_window'] in wdm_file and self.pf['hmf_model'] in wdm_file and \
+                	'_wdm_' in wdm_file:
+                    wdm_file_hmfs.append(wdm_file)
+
+            wdm_m_X_from_hmf_files = [int(hmf_file[hmf_file.find('_wdm') + 5 : hmf_file.find(\
+                '.')]) for hmf_file in wdm_file_hmfs]
+            wdm_m_X_from_hmf_files.sort()
+            #print(wdm_m_X_from_hmf_files)
+
+            closest_mass = min(wdm_m_X_from_hmf_files, key=lambda x: abs(x - m_X))
+            closest_mass_index = wdm_m_X_from_hmf_files.index(closest_mass)
+
+            if closest_mass > m_X:
+                m_X_r = closest_mass
+                m_X_l = wdm_m_X_from_hmf_files[closest_mass_index - 1]
+            elif closest_mass < m_X:
+                m_X_l = closest_mass
+                m_X_r = wdm_m_X_from_hmf_files[closest_mass_index + 1]
+            else:
+                m_X_l = int(m_X)
+                m_X_r = m_X_l + 1
+        else:
+            m_X_l = int(m_X)
+            m_X_r = m_X_l + 1
+
+        _fn = self.tab_prefix_hmf(True) + '.hdf5'
 
         if self.pf['hmf_path'] is not None:
             _path = self.pf['hmf_path'] + '/'
@@ -301,8 +336,6 @@ def _load_hmf_wdm(self): # pragma: no cover
         prefix = _fn[:_fn.find('_wdm_')]
 
         # Look for bracketing files
-        m_X_l = int(m_X)
-        m_X_r = m_X_l + 1
         fn_l = prefix + '_wdm_{:.2f}.hdf5'.format(m_X_l)
         fn_r = prefix + '_wdm_{:.2f}.hdf5'.format(m_X_r)
 
@@ -313,24 +346,39 @@ def _load_hmf_wdm(self): # pragma: no cover
         interp = True
         mass = [m_X_l, m_X_r]
         for i, fn in enumerate([fn_l, fn_r]):
+
+            # OK as long as we're not relying on interpolation
+            if not os.path.exists(_path + fn):
+                continue
+
             with h5py.File(_path + fn, 'r') as f:
 
                 tab_z = np.array(f[('tab_z')])
                 tab_M = np.array(f[('tab_M')])
                 tab_dndm = np.array(f[('tab_dndm')])
+                tab_dndm[tab_dndm==0.0] = tiny_dndm
 
                 #self.tab_k_lin = np.array(f[('tab_k_lin')])
                 #self.tab_ps_lin = np.array(f[('tab_ps_lin')])
                 #self.tab_sigma = np.array(f[('tab_sigma')])
                 #self.tab_dlnsdlnm = np.array(f[('tab_dlnsdlnm')])
 
                 tab_ngtm = np.array(f[('tab_ngtm')])
+                tab_ngtm[tab_ngtm==0.0] = tiny_dndm
                 tab_mgtm = np.array(f[('tab_mgtm')])
+                tab_mgtm[tab_mgtm==0.0] = tiny_dndm
 
                 if 'tab_MAR' in f:
-                    tab_MAR = np.array(f[('tab_MAR')])
+                    if self.pf['hmf_MAR_from_CDM']:
+                        fn_cdm = _path + prefix + '.hdf5'
+                        cdm_file = h5py.File(fn_cdm, 'r')
+                        tab_MAR = np.array(cdm_file[('tab_MAR')])
+                        cdm_file.close()
+                        print("# Loaded MAR from {}".format(fn_cdm))
+                    else:
+                        tab_MAR = np.array(f[('tab_MAR')])
                 else:
-                    print("No maR in file {}.".format(_path+fn))
+                    print("# No MAR in file {}.".format(_path+fn))
                 #self.tab_growth = np.array(f[('tab_growth')])
 
                 if m_X == mass[i]:
@@ -351,6 +399,9 @@ def _load_hmf_wdm(self): # pragma: no cover
             self.tab_mgtm = tab_mgtm
             self._tab_MAR = tab_MAR
         else:
+
+            assert len(dndm) == 2
+
             # Interpolate
             log_dndm = np.log10(dndm)
             log_ngtm = np.log10(ngtm)
@@ -366,6 +417,9 @@ def _load_hmf_wdm(self): # pragma: no cover
             self._tab_MAR = 10**(np.diff(log_tmar, axis=0).squeeze() \
                 * (m_X - m_X_l) + log_tmar[0])
 
+        if interp:
+        	print('# Finished interpolation in WDM mass dimension of HMF.')
+
     def _get_ngtm_mgtm_from_dndm(self):
 
         # Generates ngtm and mgtm with Murray's formulas
@@ -557,8 +611,10 @@ def _MF(self):
             logMmin = self.pf['hmf_logMmin']
             logMmax = self.pf['hmf_logMmax']
             dlogM = self.pf['hmf_dlogM']
-
-            from hmf.filters import SharpK, TopHat
+			#TODO FIX THIS OR REMOVE CODE
+            from hmf import filters
+            SharpK, TopHat = filters.SharpK, filters.TopHat
+            #from hmf.filters import SharpK, TopHat
             if self.pf['hmf_window'] == 'tophat':
                 # This is the default in hmf
                 window = TopHat
@@ -580,7 +636,8 @@ def _MF(self):
                 n=self.cosm.primordial_index, transfer_params=self.pars_transfer,
                 dlnk=self.pf['hmf_dlnk'], lnk_min=self.pf['hmf_lnk_min'],
                 lnk_max=self.pf['hmf_lnk_max'], hmf_params=self.pf['hmf_params'],
-                use_splined_growth=self.pf['hmf_use_splined_growth'], **xtras)
+                use_splined_growth=self.pf['hmf_use_splined_growth'],\
+                filter_params=self.pf['filter_params'], **xtras)
 
         return self._MF_
 
@@ -673,7 +730,7 @@ def TabulateHMF(self, save_MAR=True):
         MF = self._MF
 
         # Masses in hmf are really Msun / h
-        if hmf_vers < 3:
+        if hmf_vers < 3 and self.pf['hmf_wdm_mass'] is None:
             self.tab_M = self._MF.M / self.cosm.h70
         else:
             self.tab_M = self._MF.m / self.cosm.h70
@@ -1365,7 +1422,8 @@ def tab_prefix_hmf(self, with_size=False):
             s += '_{}'.format(self.pf['hmf_window'].lower())
 
         if self.pf['hmf_wdm_mass'] is not None:
-            assert self.pf['hmf_window'].lower() == 'sharpk'
+        	#TODO: For Testing, the assertion is for correct nonlinear fits.
+            #assert self.pf['hmf_window'].lower() == 'sharpk'
             s += '_wdm_{:.2f}'.format(self.pf['hmf_wdm_mass'])
 
         return s

ares/populations/Halo.py

@@ -151,4 +151,4 @@ def MGR_integrated(self, z, source=None):
 
 
 
-
+

ares/util/SetDefaultParameterValues.py

@@ -1,7 +1,7 @@
 """
 SetDefaultParameterValues.py
 
-Author: Jordan Mirocha
+Author: Jordan Mirocha / Joshua Hibbard
 Affiliation: University of Colorado at Boulder
 Created on 2010-10-19.
 
@@ -1258,6 +1258,9 @@ def HaloMassFunctionParameters():
     "hmf_wdm_mass": None,
     "hmf_wdm_interp": True,
 
+    #For various DM models
+    'hmf_dm_model': 'CDM',
+
     "hmf_cosmology_location": None,
     # PCA eigenvectors
     "hmf_pca": None,
@@ -1284,6 +1287,10 @@ def HaloMassFunctionParameters():
 
     # If a new tab_MAR should be computed when using the PCA
     "hmf_gen_MAR":False,
+
+    "filter_params" : None,
+
+    "hmf_MAR_from_CDM": True,
 
     }
 

input/hmf/generate_hmf_tables.py

@@ -6,7 +6,7 @@
 Affiliation: University of Colorado at Boulder
 Created on: Wed May  8 11:33:48 2013
 
-Description: Create lookup tables for collapsed fraction. Can be run in 
+Description: Create lookup tables for collapsed fraction. Can be run in
 parallel, e.g.,
 
     mpirun -np 4 python generate_hmf_tables.py
@@ -18,46 +18,52 @@
 
 def_kwargs = \
 {
- "hmf_model": 'PS',
+ "hmf_model": 'ST',
  "hmf_logMmin": 4,
  "hmf_logMmax": 18,
  "hmf_dlogM": 0.01,
- 
+
  "hmf_fmt": 'hdf5',
+ "hmf_table": None,
+ "hmf_wdm_mass": None,
 
- "hmf_window": 'tophat',
+ "hmf_window": 'sharpk',
 
  # Redshift sampling
- "hmf_zmin": 5.,
- "hmf_zmax": 30.,
- "hmf_dz": 0.05,
- 
- # Can do constant timestep instead of constant dz 
- #"hmf_dt": 1,
- #"hmf_tmin": 30.,
- #"hmf_tmax": 2000.,
- 
+ #"hmf_zmin": 5.,
+ #"hmf_zmax": 60.,
+ #"hmf_dz": 0.05,
+
+ # Can do constant timestep instead of constant dz
+ "hmf_dt": 1,
+ "hmf_tmin": 30.,
+ "hmf_tmax": 1000.,
+
  # Cosmology
  "cosmology_id": 'best',
  "cosmology_name": 'planck_TTTEEE_lowl_lowE',
 
+ #HMF params and filter params are for doing Aurel Schneider's 2015 paper WDM.
+ #"hmf_params" : {'a' : 1.0},	 
+ #"filter_params" : {'c' : 2.5}
+
  #"cosmology_id": 'paul',
  #"cosmology_name": 'user',
- #"sigma_8": 0.8159, 
- #'primordial_index': 0.9652, 
- #'omega_m_0': 0.315579, 
- #'omega_b_0': 0.0491, 
+ #"sigma_8": 0.8159,
+ #'primordial_index': 0.9652,
+ #'omega_m_0': 0.315579,
+ #'omega_b_0': 0.0491,
  #'hubble_0': 0.6726,
- #'omega_l_0': 1. - 0.315579, 
- 
+ #'omega_l_0': 1. - 0.315579,
+
 }
 
 ##
 
 kwargs = def_kwargs.copy()
 kwargs.update(ares.util.get_cmd_line_kwargs(sys.argv))
 
-hmf = ares.physics.HaloMassFunction(hmf_analytic=False, 
+hmf = ares.physics.HaloMassFunction(hmf_analytic=False,
     hmf_load=False, **kwargs)
 
 hmf.info
@@ -66,6 +72,3 @@
     hmf.SaveHMF(fmt=kwargs['hmf_fmt'], clobber=False)
 except IOError as err:
     print(err)
-
-
-