From 47ae5c6fc4187dcbd9fc27a65110d2087a709354 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Mon, 14 May 2018 02:02:02 -0500 Subject: [PATCH 1/6] Modernizing IMF formulation; adding factory call; adding testing function --- tests/test_imf.py | 43 ++++++++ ugali/analysis/imf.py | 168 ++++++++++++++++++------------ ugali/config/config_dr10_gal.yaml | 2 +- ugali/config/config_sva1.yaml | 2 +- ugali/isochrone/model.py | 4 +- 5 files changed, 149 insertions(+), 70 deletions(-) create mode 100644 tests/test_imf.py diff --git a/tests/test_imf.py b/tests/test_imf.py new file mode 100644 index 0000000..e75ee87 --- /dev/null +++ b/tests/test_imf.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python +""" +Generic python script. +""" +__author__ = "Alex Drlica-Wagner" +import numpy as np + +import ugali.analysis.imf +import ugali.isochrone + +MASSES = np.array([0.1, 0.5, 1.0, 10.0]) +SEED = 1 + +def test_imf(): + dn_dm = ugali.analysis.imf.chabrierIMF(MASSES) + np.testing.assert_allclose(dn_dm, + [1.29919665, 0.66922584, 0.3666024, 0.01837364], + rtol=1e-6) + + #imf = ugali.analysis.imf.IMF('chabrier') + imf = ugali.analysis.imf.factory('Chabrier2003') + masses = imf.sample(3, seed=SEED) + np.testing.assert_allclose(masses,[0.22122041, 0.4878909, 0.10002029],rtol=1e-6) + + imf = ugali.analysis.imf.Chabrier2003() + masses = imf.sample(3, seed=SEED) + np.testing.assert_allclose(masses,[0.22122041, 0.4878909, 0.10002029],rtol=1e-6) + + iso = ugali.isochrone.Bressan2012(imf_type='Chabrier2003') + masses = iso.imf.sample(3,seed=SEED) + np.testing.assert_allclose(masses,[0.22122041, 0.4878909, 0.10002029],rtol=1e-6) + + integral = iso.imf.integrate(0.1,2.0) + np.testing.assert_allclose(integral,0.94961632593) + + import pdb; pdb.set_trace() + +if __name__ == "__main__": + import argparse + parser = argparse.ArgumentParser(description=__doc__) + args = parser.parse_args() + + test_imf() diff --git a/ugali/analysis/imf.py b/ugali/analysis/imf.py index 2f3184e..1e179d6 100644 --- a/ugali/analysis/imf.py +++ b/ugali/analysis/imf.py @@ -1,8 +1,7 @@ """ Classes to handle initial mass functions (IMFs). """ - -import numpy +from abc import abstractmethod import numpy as np import scipy.interpolate @@ -13,35 +12,25 @@ ############################################################ -class IMF: +class IMF(object): """ Base class for initial mass functions (IMFs). """ - def __init__(self, type='chabrier'): - """ - Initialize an instance of an initial mass function. - """ - - self.type = type - - if self.type == 'chabrier': - self.pdf = chabrierIMF - else: - logger.warn('initial mass function type %s not recognized'%(self.type)) - - def integrate(self, mass_min, mass_max, log_mode=True, weight=False, steps=10000): - """ - Numerically integrate initial mass function. + def integrate(self, mass_min, mass_max, log_mode=True, weight=False, steps=1e4): + """ Numerically integrate the IMF. - INPUTS: - mass_min: minimum mass bound for integration (solar masses) - mass_max: maximum mass bound for integration (solar masses) - log_mode[True]: use logarithmic steps in stellar mass as oppose to linear + Parameters: + ----------- + mass_min: minimum mass bound for integration (solar masses) + mass_max: maximum mass bound for integration (solar masses) + log_mode[True]: use logarithmic steps in stellar mass as oppose to linear weight[False]: weight the integral by stellar mass steps: number of numerical integration steps - OUTPUT: - result of integral + + Returns: + -------- + result of integral """ if log_mode: d_log_mass = (np.log10(mass_max) - np.log10(mass_min)) / float(steps) @@ -63,59 +52,106 @@ def integrate(self, mass_min, mass_max, log_mode=True, weight=False, steps=10000 def sample(self, n, mass_min=0.1, mass_max=10., steps=10000, seed=None): """ - Sample n initial mass values between mass_min and mass_max, following the IMF distribution. + Sample initial mass values between mass_min and mass_max, + following the IMF distribution. + + ADW: Should this be `sample` or `simulate`? + + Parameters: + ----------- + n : number of samples to draw + mass_min : minimum mass to sample from + mass_max : maximum mass to sample from + steps : number of steps for isochrone sampling + seed : random seed (passed to np.random.seed) + + Returns: + -------- + mass : array of randomly sampled mass values """ if seed is not None: np.random.seed(seed) d_mass = (mass_max - mass_min) / float(steps) - mass = numpy.linspace(mass_min, mass_max, steps) - cdf = numpy.insert(numpy.cumsum(d_mass * self.pdf(mass[1:], log_mode=False)), 0, 0.) + mass = np.linspace(mass_min, mass_max, steps) + cdf = np.insert(np.cumsum(d_mass * self.pdf(mass[1:], log_mode=False)), 0, 0.) cdf = cdf / cdf[-1] f = scipy.interpolate.interp1d(cdf, mass) - return f(numpy.random.uniform(size=n)) + return f(np.random.uniform(size=n)) -############################################################ - -def chabrierIMF(mass, log_mode=True, a=1.31357499301): - """ - Chabrier initial mass function. + @abstractmethod + def pdf(cls): pass + +class Chabrier2003(IMF): + """ Initial mass function from Chabrier (2003): "Galactic Stellar and Substellar Initial Mass Function" Chabrier PASP 115:763-796 (2003) https://arxiv.org/abs/astro-ph/0304382 + """ + + @classmethod + def pdf(cls, mass, log_mode=True, a=1.31357499301): + """ PDF for the Chabrier IMF. + + The functional form and coefficients are described in Eq 17 + and Tab 1 of Chabrier (2003): + + m <= 1 Msun: E(log m) = A1*exp(-(log m - log m_c)^2 / 2 sigma^2) + m > 1 Msun: E(log m) = A2 * m^-x + + A1 = 1.58 : normalization [ log(Msun)^-1 pc^-3] + m_c = 0.079 [Msun] + sigma = 0.69 + A2 = 4.43e-2 + x = 1.3 + + We redefine a = A1, A2 = a * b; + + Chabrier set's his normalization + + Parameters: + ----------- + mass: stellar mass (solar masses) + log_mode[True]: return number per logarithmic mass range, i.e., dN/dlog(M) + a[1.31357499301]: normalization; normalized by default to the mass interval 0.1--100 solar masses + + Returns: + -------- + number per (linear or log) mass bin, i.e., dN/dM or dN/dlog(M) where mass unit is solar masses + """ + log_mass = np.log10(mass) + + # Constants from Chabrier 2003 + m_c = 0.079 + sigma = 0.69 + x = 1.3 + + # This value is required so that the two components match at 1 Msun + b = 0.279087531047 + + dn_dlogm = ((log_mass <= 0) * a * np.exp(-(log_mass - np.log10(m_c))**2 / (2 * (sigma**2)))) + ((log_mass > 0) * a * b * mass**(-x)) + + if log_mode: + # Number per logarithmic mass range, i.e., dN/dlog(M) + return dn_dlogm + else: + # Number per linear mass range, i.e., dN/dM + return dn_dlogm / (mass * np.log(10)) + + +class Kroupa2002(IMF): pass +class Salpeter2002(IMF): pass + +def factory(name, **kwargs): + from ugali.utils.factory import factory + return factory(name, module=__name__, **kwargs) + +kernelFactory = factory - The form and coefficients are described in Equation 17 and Table 1: - m <= 1 Msun: E(log m) = A1*exp(-(log m - log m_c)^2 / 2 sigma^2) - m > 1 Msun: E(log m) = A2 * m^-x - - A1 = 1.58 : normalization [ log(Msun)^-1 pc^-3] - m_c = 0.079 [Msun] - sigma = 0.69 - A2 = 4.43e-2 - x = 1.3 - - We redefine a = A1, A2 = a * b; - - Chabrier set's his normalization - - Parameters: - ----------- - mass: stellar mass (solar masses) - log_mode[True]: return number per logarithmic mass range, i.e., dN/dlog(M) - a[1.31357499301]: normalization; normalized by default to the mass interval 0.1--100 solar masses - Returns: - -------- - number per (linear or log) mass bin, i.e., dN/dM or dN/dlog(M) where mass unit is solar masses - """ - log_mass = np.log10(mass) - b = 0.279087531047 # This value is required so that the two components match at 1 Msun - if log_mode: - # Number per logarithmic mass range, i.e., dN/dlog(M) - return ((log_mass <= 0) * a * np.exp(-(log_mass - np.log10(0.079))**2 / (2 * (0.69**2)))) + \ - ((log_mass > 0) * a * b * mass**(-1.3)) - else: - # Number per linear mass range, i.e., dN/dM - return (((log_mass <= 0) * a * np.exp(-(log_mass - np.log10(0.079))**2 / (2 * (0.69**2)))) + \ - ((log_mass > 0) * a * b * mass**(-1.3))) / (mass * np.log(10)) +############################################################ + +def chabrierIMF(mass, log_mode=True, a=1.31357499301): + """ Backward compatible wrapper around Chabrier2003.pdf """ + return Chabrier2003.pdf(mass,log_mode=log_mode,a=a) ############################################################ diff --git a/ugali/config/config_dr10_gal.yaml b/ugali/config/config_dr10_gal.yaml index 71bd362..c098306 100644 --- a/ugali/config/config_dr10_gal.yaml +++ b/ugali/config/config_dr10_gal.yaml @@ -179,7 +179,7 @@ isochrone: mag_1_field: g mag_2_field: r stage_field: stage - imf: chabrier + imf: Chabrier2003 horizontal_branch_stage : 'BHeb' #horizontal_branch_dispersion: null #horizontal_branch_dispersion: [0.0] # mag diff --git a/ugali/config/config_sva1.yaml b/ugali/config/config_sva1.yaml index f43e210..130286c 100644 --- a/ugali/config/config_sva1.yaml +++ b/ugali/config/config_sva1.yaml @@ -172,7 +172,7 @@ isochrone: mag_1_field: g mag_2_field: r stage_field: stage - imf: chabrier + imf: Chabrier2003 horizontal_branch_stage : 'BHeb' #horizontal_branch_dispersion: null #horizontal_branch_dispersion: [0.0] # mag diff --git a/ugali/isochrone/model.py b/ugali/isochrone/model.py index 36e3cd0..cdec43e 100644 --- a/ugali/isochrone/model.py +++ b/ugali/isochrone/model.py @@ -84,7 +84,7 @@ class IsochroneModel(Model): ('band_1','g','Field name for magnitude one'), ('band_2','r','Field name for magnitude two'), ('band_1_detection',True,'Band one is detection band'), - ('imf_type','chabrier','Initial mass function'), + ('imf_type','Chabrier2003','Initial mass function'), ('hb_stage',None,'Horizontal branch stage name'), ('hb_spread',0.0,'Intrinisic spread added to horizontal branch'), ) @@ -100,7 +100,7 @@ def _setup(self, **kwargs): for k,v in list(defaults.items()): setattr(self,k,v) - self.imf = ugali.analysis.imf.IMF(defaults['imf_type']) + self.imf = ugali.analysis.imf.factory(defaults['imf_type']) self.index = None def _parse(self,filename): From f18ecc0a0246fd5a5b372fbd27c87ab150ad5909 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Mon, 14 May 2018 02:14:21 -0500 Subject: [PATCH 2/6] remove pdb --- tests/test_imf.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_imf.py b/tests/test_imf.py index e75ee87..9c5daf9 100644 --- a/tests/test_imf.py +++ b/tests/test_imf.py @@ -33,8 +33,6 @@ def test_imf(): integral = iso.imf.integrate(0.1,2.0) np.testing.assert_allclose(integral,0.94961632593) - import pdb; pdb.set_trace() - if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description=__doc__) From ac167f6c3d896ebee70b3d07de6e5aa0c24cd8e0 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Mon, 14 May 2018 19:39:38 -0500 Subject: [PATCH 3/6] Adding Salpeter1955 and Kroupa2001 IMFs --- tests/test_imf.py | 12 +++++- ugali/analysis/imf.py | 99 +++++++++++++++++++++++++++++++++++-------- 2 files changed, 92 insertions(+), 19 deletions(-) diff --git a/tests/test_imf.py b/tests/test_imf.py index 9c5daf9..fcc2529 100644 --- a/tests/test_imf.py +++ b/tests/test_imf.py @@ -10,6 +10,7 @@ MASSES = np.array([0.1, 0.5, 1.0, 10.0]) SEED = 1 +IMFS = ['Salpeter1955','Kroupa2001','Chabrier2003'] def test_imf(): dn_dm = ugali.analysis.imf.chabrierIMF(MASSES) @@ -33,9 +34,18 @@ def test_imf(): integral = iso.imf.integrate(0.1,2.0) np.testing.assert_allclose(integral,0.94961632593) +def test_imf_norm(): + # All of the IMFs have normalized integrals from 0.1 to 100 Msun + for n in IMFS: + print n + imf = ugali.analysis.imf.factory(name=n) + np.testing.assert_allclose(imf.integrate(0.1,100,steps=int(1e4)), + 1.0,rtol=1e-3) + if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description=__doc__) args = parser.parse_args() - test_imf() + #test_imf() + test_imf_norm() diff --git a/ugali/analysis/imf.py b/ugali/analysis/imf.py index 1e179d6..e585538 100644 --- a/ugali/analysis/imf.py +++ b/ugali/analysis/imf.py @@ -1,5 +1,7 @@ """ Classes to handle initial mass functions (IMFs). + +https://github.com/keflavich/imf """ from abc import abstractmethod import numpy as np @@ -7,9 +9,6 @@ from ugali.utils.logger import logger -# ADW: TODO - This needs to be modernized -# Also add Kroupa and Salpeter IMFs... - ############################################################ class IMF(object): @@ -17,8 +16,12 @@ class IMF(object): Base class for initial mass functions (IMFs). """ + def __call__(self, mass, **kwargs): + """ Call the pdf of the mass function """ + return self.pdf(mass,**kwargs) + def integrate(self, mass_min, mass_max, log_mode=True, weight=False, steps=1e4): - """ Numerically integrate the IMF. + """ Numerical Riemannn integral of the IMF (stupid simple). Parameters: ----------- @@ -78,7 +81,7 @@ def sample(self, n, mass_min=0.1, mass_max=10., steps=10000, seed=None): return f(np.random.uniform(size=n)) @abstractmethod - def pdf(cls): pass + def pdf(cls, mass, **kwargs): pass class Chabrier2003(IMF): @@ -89,7 +92,7 @@ class Chabrier2003(IMF): """ @classmethod - def pdf(cls, mass, log_mode=True, a=1.31357499301): + def pdf(cls, mass, log_mode=True): """ PDF for the Chabrier IMF. The functional form and coefficients are described in Eq 17 @@ -105,14 +108,14 @@ def pdf(cls, mass, log_mode=True, a=1.31357499301): x = 1.3 We redefine a = A1, A2 = a * b; - - Chabrier set's his normalization - + + The normalization is set so that the IMF integrates to 1 over + the mass range from 0.1 Msun to 100 Msun + Parameters: ----------- mass: stellar mass (solar masses) log_mode[True]: return number per logarithmic mass range, i.e., dN/dlog(M) - a[1.31357499301]: normalization; normalized by default to the mass interval 0.1--100 solar masses Returns: -------- @@ -123,31 +126,91 @@ def pdf(cls, mass, log_mode=True, a=1.31357499301): # Constants from Chabrier 2003 m_c = 0.079 sigma = 0.69 - x = 1.3 + x = 1.3 + # This value is set to normalize integral from 0.1 to 100 Msun + a=1.31357499301 # This value is required so that the two components match at 1 Msun b = 0.279087531047 - dn_dlogm = ((log_mass <= 0) * a * np.exp(-(log_mass - np.log10(m_c))**2 / (2 * (sigma**2)))) + ((log_mass > 0) * a * b * mass**(-x)) - + dn_dlogm = ((log_mass <= 0) * a * np.exp(-(log_mass - np.log10(m_c))**2 / (2 * (sigma**2)))) + dn_dlogm += ((log_mass > 0) * a * b * mass**(-x)) + if log_mode: # Number per logarithmic mass range, i.e., dN/dlog(M) return dn_dlogm else: # Number per linear mass range, i.e., dN/dM return dn_dlogm / (mass * np.log(10)) - -class Kroupa2002(IMF): pass -class Salpeter2002(IMF): pass +class Kroupa2001(IMF): + """ IMF from Kroupa (2001): + + "On the variation of the initial mass function" + MNRAS 322:231 (2001) + https://arxiv.org/abs/astro-ph/0009005 + """ + + @classmethod + def pdf(cls, mass, log_mode=True): + """ PDF for the Kroupa IMF. + + Normalization is set over the mass range from 0.1 Msun to 100 Msun + """ + log_mass = np.log10(mass) + # From Eq 2 + mb = mbreak = [0.08, 0.5] # Msun + a = alpha = [0.3, 1.3, 2.3] # alpha + + # Normalization set from 0.1 -- 100 Msun + norm = 0.27947743949440446 + + b = 1./norm + c = b * mbreak[0]**(alpha[1]-alpha[0]) + d = c * mbreak[1]**(alpha[2]-alpha[1]) + + dn_dm = b * (mass < 0.08) * mass**(-alpha[0]) + dn_dm += c * (0.08 <= mass) * (mass < 0.5) * mass**(-alpha[1]) + dn_dm += d * (0.5 <= mass) * mass**(-alpha[2]) + if log_mode: + # Number per logarithmic mass range, i.e., dN/dlog(M) + return dn_dm * (mass * np.log(10)) + else: + # Number per linear mass range, i.e., dN/dM + return dn_dm + +class Salpeter1955(IMF): + """ IMF from Salpeter (1955): + "The Luminosity Function and Stellar Evolution" + ApJ 121, 161S (1955) + http://adsabs.harvard.edu/abs/1955ApJ...121..161S + """ + + @classmethod + def pdf(cls, mass, log_mode=True): + """ PDF for the Salpeter IMF. + + Value of 'a' is set to normalize the IMF to 1 between 0.1 and 100 Msun + """ + alpha = 2.35 + + a = 0.060285569480482866 + dn_dm = a * mass**(-alpha) + + if log_mode: + # Number per logarithmic mass range, i.e., dN/dlog(M) + return dn_dm * (mass * np.log(10)) + else: + # Number per linear mass range, i.e., dN/dM + return dn_dm + def factory(name, **kwargs): from ugali.utils.factory import factory return factory(name, module=__name__, **kwargs) -kernelFactory = factory +imfFactory = factory - ############################################################ def chabrierIMF(mass, log_mode=True, a=1.31357499301): From c6e088ece0b07ce75703f738206af96ebb8f41f7 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 15 May 2018 00:28:59 -0500 Subject: [PATCH 4/6] removing 'a' parameter --- ugali/analysis/imf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ugali/analysis/imf.py b/ugali/analysis/imf.py index e585538..81c391b 100644 --- a/ugali/analysis/imf.py +++ b/ugali/analysis/imf.py @@ -213,8 +213,8 @@ def factory(name, **kwargs): ############################################################ -def chabrierIMF(mass, log_mode=True, a=1.31357499301): +def chabrierIMF(mass, log_mode=True): """ Backward compatible wrapper around Chabrier2003.pdf """ - return Chabrier2003.pdf(mass,log_mode=log_mode,a=a) + return Chabrier2003.pdf(mass,log_mode=log_mode) ############################################################ From 5be681956cd0d8653b918138221237e5761c2b27 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 15 May 2018 01:07:04 -0500 Subject: [PATCH 5/6] catching download errors --- tests/test_isochrone.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_isochrone.py b/tests/test_isochrone.py index 679ccc0..5204938 100644 --- a/tests/test_isochrone.py +++ b/tests/test_isochrone.py @@ -142,9 +142,9 @@ def test_download(): try: iso.download(outdir='./tmp/'+name.lower(),force=True) except URLError as e: - print("Server is down") - print(e) - + logger.error("Server is down:\n%s"%str(e)) + except RuntimeError as e: + logger.error(str(e)) if __name__ == "__main__": import argparse From 05c2f7cfd6c9f07f5b18f3d76488ae6c6b66421b Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 15 May 2018 01:27:45 -0500 Subject: [PATCH 6/6] py3 print --- tests/test_imf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_imf.py b/tests/test_imf.py index fcc2529..f5850e0 100644 --- a/tests/test_imf.py +++ b/tests/test_imf.py @@ -37,7 +37,7 @@ def test_imf(): def test_imf_norm(): # All of the IMFs have normalized integrals from 0.1 to 100 Msun for n in IMFS: - print n + print(n) imf = ugali.analysis.imf.factory(name=n) np.testing.assert_allclose(imf.integrate(0.1,100,steps=int(1e4)), 1.0,rtol=1e-3)