Error when trying to make_cluster with mmin > break1 in Kroupa mass function

[hcferguson]

Thanks for this very useful package!

In simulating nearby dwarf galaxies, I'd like to save time in make_cluster by setting a lower-mass limit that's pretty high. This worked with an earlier version of this package (or at least didn't throw an error). Now I'm getting an error having to do with checking where the breaks in the power law distribution are.

It works if I do cl = imf.make_cluster(1000,massfunc='kroupa') but throws an error if I do cl = imf.make_cluster(1000,massfunc='kroupa',mmin=1.0).

The weird thing is that if I then try to repeat the first command again cl = imf.make_cluster(1000,massfunc='kroupa'), it now throws a different error.

Traceback from the first error:

cl = imf.make_cluster(1000,massfunc='kroupa',mmin=1.0)
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-3-878462dd2f07> in <module>
----> 1 cl = imf.make_cluster(1000,massfunc='kroupa',mmin=1.0)

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/imf.py in make_cluster(mcluster, massfunc, verbose, silent, tolerance, stop_criterion, mmax, mmin)
    451     if mmin is not None and hasattr(mfc, 'mmin') and mfc.mmin != mmin:
    452         orig_mmin = mfc.mmin
--> 453         mfc.mmin = mmin
    454     if mmax is not None and hasattr(mfc, 'mmax') and mfc.mmax != mmax:
    455         orig_mmax = mfc.mmax

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/imf.py in mmin(self, value)
    129     @mmin.setter
    130     def mmin(self, value):
--> 131         self.distr.m1 = value
    132 
    133     @property

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/distributions.py in m1(self, value)
    152     def m1(self, value):
    153         self.breaks[0] = value
--> 154         self._calcpows()
    155         self._calcweights()
    156 

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/distributions.py in _calcpows(self)
    169         pows = []
    170         for ii in range(nsegm):
--> 171             pows.append(PowerLaw(self.slopes[ii], self.breaks[ii], self.breaks[ii + 1]))
    172         self.pows = pows
    173         self.nsegm = nsegm

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/distributions.py in __init__(self, slope, m1, m2)
     88         self.m1 = float(m1)
     89         self.m2 = float(m2)
---> 90         assert (m1 < m2)
     91         assert (m1 > 0)
     92         assert (m1 != -1)

AssertionError: 

Traceback from the second time

cl = imf.make_cluster(1000,massfunc='kroupa')
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-4-cac326b32cec> in <module>
----> 1 cl = imf.make_cluster(1000,massfunc='kroupa')

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/imf.py in make_cluster(mcluster, massfunc, verbose, silent, tolerance, stop_criterion, mmax, mmin)
    460         assert expected_mass > 0
    461     else:
--> 462         expected_mass = mfc.m_integrate(mfc.mmin, mmax)[0]
    463         assert expected_mass > 0
    464         expectedmass_cache[(massfunc, mfc.mmin, mmax)] = expected_mass

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/imf.py in m_integrate(self, mlow, mhigh, numerical, **kwargs)
    186                                                    -self.p3+1],
    187                                                   [self.mmin, self.break1,
--> 188                                                    self.break2, self.mmax])
    189             ratio = distr1.pdf(self.break1)/self.distr.pdf(self.break1)/self.break1
    190             return ((distr1.cdf(mhigh)-distr1.cdf(mlow))/ratio, 0)

~/anaconda3/envs/astroconda37/lib/python3.7/site-packages/imf-0.1.dev151+g9a33b9e-py3.7.egg/imf/distributions.py in __init__(self, slopes, breaks)
    139         """
    140         assert (len(slopes) == len(breaks) - 1)
--> 141         assert ((np.diff(breaks) > 0).all())
    142         self.slopes = slopes
    143         self.breaks = breaks

AssertionError: 

[hcferguson]

I think I see why the weirdness of the second error is happening. When the first attempt fails, make_cluster never makes it down to the lines where it restores the original values of mmin and mmax:

    if 'orig_mmin' in locals():
        mfc.mmin = orig_mmin
    if 'orig_mmax' in locals():
        mfc.mmax = orig_mmax

[keflavich]

Oof, that's an ugly error. Can I suggest, though, that in this case you just want a Salpeter (powerlaw) IMF instead of a Kroupa, since you're only using one powerlaw?

[keflavich]

This error still occurs; if you want to use a Kroupa MF with a minimum mass above the first break, you have to define it custom following the example in the code:

class Kroupa(MassFunction):
    # kroupa = BrokenPowerLaw(breaks={0.08: -0.3, 0.5: 1.3, 'last': 2.3}, mmin=0.03, mmax=120)
    default_mmin = 0.03
    default_mmax = 120

    def __init__(self,
                 mmin=default_mmin,
                 mmax=default_mmax,
                 p1=0.3,
                 p2=1.3,
                 p3=2.3,
                 break1=0.08,
                 break2=0.5):
        """
        The Kroupa IMF with two power-law breaks, p1 and p2. See __call__ for
        details.
        """
        super().__init__(mmin=mmin, mmax=mmax)

        self.p1 = p1
        self.p2 = p2
        self.p3 = p3
        self.break1 = break1
        self.break2 = break2
        self.distr = distributions.BrokenPowerLaw([-p1, -p2, -p3],
                                                  [self.mmin, break1, break2, self.mmax])
        self.normfactor = 1