Mass function

massFunc.py

def mass_function(data, logged=False,
                  nbins=35, verbose=True):
    """
    Calculates a mass function from data, using 3 methods.

    Parameters
    ----------
    data: array

    logged: boolean
        if False, will take log inside function

    nbins: int

    verbose: boolean


    Return
    ------
    mf / dm:
        dN/dlnM, differential mass function
    mbin:
        mass bins

    """

    nbins = int(nbins)

    #if log have been taken from the data or not
    if not logged:
        dat = np.log(data)
    else:
        dat = data
    del data

    #number of object
    nobj = len(dat)

    #bins
    mmax = dat.max()
    mmin = dat.min()
    dm = (mmax - mmin) / float(nbins)
    #mid points
    mbin = (np.arange(nbins) + 0.5) * dm + mmin

    # count up mass in each bin
    mf = np.zeros(nbins)
    # which bin each data point belongs to
    ibin = np.floor((dat - mmin) / dm)
    # make a mask to keep bins given user-defind nbins
    mask = (ibin >= 0) & (ibin < nbins)
    # sum up each bin
    wght = np.ones(len(ibin[mask]))
    for i in range(nbins):
        mf[i] = np.sum(wght[np.array(ibin[mask], dtype=int) == i])
    mf = np.array(mf)

    if not logged:
        mbin = np.e ** mbin

    if verbose:
        print 'Number of object = %i' % nobj
        print 'dlnM =', dm
        print 'min = %f, max = %f' % (mmin, mmax)
        print 'Results:\n', mbin, mf / dm

    # alternatively, we can calculating using np.digitize....
    id_bin = np.digitize(dat, mbin)
    mf = np.bincount(id_bin)
    print "here: ", mf/dm

    # calculate mass functions, using np.histogram:
    # counts_per_bin, bin_edges
    mf, edges = np.histogram(dat, nbins,
                            range=(mmin, mmax))

    mbin = (edges[1:] + edges[:-1]) / 2.
    dm = edges[1] - edges[0]

    if not logged:
        mbin = np.e ** mbin

    if verbose:
        print 'From np.histogram: '
        print 'Number of object = %i' % nobj
        print 'dlnM =', dm
        print 'min = %f, max = %f' % (mmin, mmax)
        print 'Results:\n', mbin, mf / dm

    return mbin, mf / dm