Finish 3.0.3

# Conflicts:
#	CHANGELOG.rst

CHANGELOG.rst

@@ -8,8 +8,18 @@ dev
 - Added `get_dependencies` method to _Framework, to enable finding all parameters that a quantity depends on.
 
 
-v3.0.2 [3rd Nov]
-----------------
+v3.0.3 [1st Dec 2017]
+---------------------
+**Bugfixes**
+
+- Fixed usage of deprecated MsolMass in wdm
+- Fixed Bhattachrya fitting function (thanks to Benedikt Diemer!)
+- Fixed typo in Watson fitting function (thanks to Benedikt Diemer!)
+- Update cosmo test to use new Astropy constants.
+- Fixed issue with sampling function where zeros in ngtm would yield an error.
+
+v3.0.2 [3rd Nov 2017]
+---------------------
 **Bugfixes**
 
 - Changed parameter checks on instantiation to post-conversion.

hmf/cosmo.py

@@ -101,7 +101,7 @@ def mean_density0(self):
         """
         Mean density of universe at z=0, [Msun h^2 / Mpc**3]
         """
-        return (self.cosmo.Om0 * self.cosmo.critical_density0 / self.cosmo.h ** 2).to(u.solMass/u.Mpc**3).value
+        return (self.cosmo.Om0 * self.cosmo.critical_density0 / self.cosmo.h ** 2).to(u.Msun/u.Mpc**3).value
 
 def get_cosmo(name):
     """

hmf/fitting_functions.py

@@ -558,7 +558,7 @@ def cutmask(self):
 class Angulo(FittingFunction):
     req_mass = True
     _ref = """Angulo, R. E., et al., 2012. arXiv:1203.3216v1"""
-    _eq = r"$A \left(\frac{d}{\sigma}\right)^b \exp(-c/\sigma^2)$"
+    _eq = r"$A \left[\left(\frac{d}{\sigma}\right)^b + 1 \right] \exp(-c/\sigma^2)$"
     __doc__ = _makedoc(FittingFunction._pdocs, "Angulo", "Ang", _eq, _ref)
     _defaults = {"A": 0.201, "b": 1.7, "c": 1.172, "d": 2.08}
 
@@ -587,7 +587,7 @@ def fsigma(self):
         c = self.params['c']
         d = self.params['d']
 
-        return 0.201 * ((d/self.sigma)**b+1) * np.exp(-c/self.sigma**2)
+        return A * ((d/self.sigma)**b+1) * np.exp(-c/self.sigma**2)
 
     @property
     def cutmask(self):
@@ -644,7 +644,7 @@ class Watson(FittingFunction):
     sim_definition.halo_overdensity = 178.0
 
     _defaults = {"C_a": 0.023, "d_a": 0.456, "d_b": 0.139, "p": 0.072, "q": 2.13,
-                 "A_0": 0.194, "alpha_0": 2.267, "beta_0": 1.805, "gamma_0": 1.287,
+                 "A_0": 0.194, "alpha_0": 1.805, "beta_0": 2.267, "gamma_0": 1.287,
                  "z_hi": 6, "A_hi": 0.563, "alpha_hi": 0.874, "beta_hi": 3.810, "gamma_hi": 1.453,
                  "A_a": 1.097, "A_b": 3.216, "A_c": 0.074,
                  "alpha_a": 3.136, "alpha_b": 3.058, "alpha_c": 2.349,
@@ -765,7 +765,7 @@ class Bhattacharya(SMT):
     req_z = True
     req_mass = True
 
-    _eq = r"f_{\rm SMT}(\sigma) (\nu\sqrt{a})^q"
+    _eq = r"f_{\rm SMT}(\sigma) (\nu\sqrt{a})^{q-1}"
     _ref = """Bhattacharya, S., et al., May 2011. ApJ 732 (2), 122. http://labs.adsabs.harvard.edu/ui/abs/2011ApJ...732..122B"""
     __doc__ = _makedoc(FittingFunction._pdocs, "Bhattacharya", "Btc", _eq, _ref)
     _defaults = {"A_a": 0.333, "A_b": 0.11, "a_a": 0.788, "a_b": 0.01, "p": 0.807, "q": 1.795}
@@ -810,7 +810,7 @@ def fsigma(self):
             The function :math:`f(\sigma)\equiv\nu f(\nu)` defined on ``pert.M``
         """
         vfv = super(Bhattacharya, self).fsigma
-        return vfv*self.params['a']**(self.params['q']/2)*self.nu**(self.params['q'] - 1)
+        return vfv*(np.sqrt(self.params['a']) * self.nu)**(self.params['q'] - 1)
 
     @property
     def cutmask(self):

hmf/sample.py

@@ -9,7 +9,8 @@
 
 def _prepare_mf(log_mmin, **mf_kwargs):
     h = hmf.MassFunction(Mmin=log_mmin, **mf_kwargs)
-    icdf = _spline((h.ngtm / h.ngtm[0])[::-1], np.log10(h.m[::-1]), k=3)
+    mask = h.ngtm>0
+    icdf = _spline((h.ngtm[mask] / h.ngtm[0])[::-1], np.log10(h.m[mask][::-1]), k=3)
 
     return icdf, h
 
@@ -51,12 +52,6 @@ def sample_mf(N,log_mmin,
     hmf : `hmf.MassFunction` instance
         The instance used to define the mass function.
 
-    Notes
-    -----
-    `Mmax` is a free parameter to be sent to `MassFunction`. However, if set
-    at or above 18, the routine will fail. It should be set so that dndm is
-    very small at Mmax, but not numerically 0.
-
     Examples
     --------
     Simplest example:

hmf/wdm.py

@@ -47,7 +47,7 @@ class attribute.
     def __init__(self, mx, cosmo = Planck15, z =0, **model_params):
         self.mx = mx
         self.cosmo = cosmo
-        self.rho_mean = (1+z)**3 * (self.cosmo.Om0 * self.cosmo.critical_density0 / self.cosmo.h ** 2).to(u.MsolMass / u.Mpc ** 3).value * 1e6
+        self.rho_mean = (1+z)**3 * (self.cosmo.Om0 * self.cosmo.critical_density0 / self.cosmo.h ** 2).to(u.solMass / u.Mpc ** 3).value * 1e6
         self.Oc0 = cosmo.Om0 - cosmo.Ob0
 
         super(WDM, self).__init__(**model_params)

tests/test_cosmo.py

@@ -8,6 +8,7 @@
 from hmf.cosmo import Cosmology, WMAP7
 import numpy as np
 
+
 def eq(actual, expected):
     return abs(actual - expected) < 0.000001
 
@@ -23,10 +24,12 @@ def eq(actual, expected):
 #     for bit in bits:
 #         assert bit in c.pycamb_dict
 
+
 def test_string_cosmo():
     c = Cosmology(cosmo_model="WMAP7")
     assert c.cosmo.Ob0 > 0
 
+
 class TestUpdate():
     def __init__(self):
         self.c = Cosmology(cosmo_model="Planck13")
@@ -35,8 +38,8 @@ def test_cosmo_model(self):
         self.c.update(cosmo_model=WMAP7)
 
         assert self.c.cosmo.Om0 == 0.272
-        print((self.c.mean_density0))
-        assert np.isclose(self.c.mean_density0, 75468972351.60081)
+        print(self.c.mean_density0)
+        assert np.isclose(self.c.mean_density0, 75489962610.27452, atol=1e-3) # this number *can* change when updated constants are used.
 
     def test_cosmo_params(self):
         self.c.update(cosmo_params={"H0":0.6})

tests/test_sample.py

@@ -23,6 +23,7 @@ def test_mmax_big():
     #raises ValueError because ngtm=0 exactly at m=18
     #due to hard limit of integration in integrate_hmf.
     np.random.seed(12345)
+
     m,h = sample_mf(1e5,11,transfer_model="EH",Mmax=18)
     #centres,hist = dndm_from_sample(m,1e5/h.ngtm[0])
     #print centres,hist