utils

photXclus.utils

cat2hpx(lon, lat, nside, radec=True)

Convert a catalogue to a HEALPix map of number counts per resolution element.

Parameters:

Name	Type	Description	Default
lon	(ndarray, ndarray)	Coordinates of the sources in degree. If radec=True, assume input is in the icrs coordinate system. Otherwise assume input is glon, glat	required
lat	(ndarray, ndarray)	Coordinates of the sources in degree. If radec=True, assume input is in the icrs coordinate system. Otherwise assume input is glon, glat	required
nside	int	HEALPix nside of the target map	required
radec	bool	Switch between R.A./Dec and glon/glat as input coordinate system.	True

Return

hpx_map : ndarray HEALPix map of the catalogue number counts in Galactic coordinates

Source code in photXclus/utils.py

def cat2hpx(lon, lat, nside, radec=True):
    """
    Convert a catalogue to a HEALPix map of number counts per resolution
    element.

    Parameters
    ----------
    lon, lat : (ndarray, ndarray)
        Coordinates of the sources in degree. If radec=True, assume input is in the icrs
        coordinate system. Otherwise assume input is glon, glat

    nside : int
        HEALPix nside of the target map

    radec : bool
        Switch between R.A./Dec and glon/glat as input coordinate system.

    Return
    ------
    hpx_map : ndarray
        HEALPix map of the catalogue number counts in Galactic coordinates

    """

    npix = hp.nside2npix(nside)

    if radec:
        eq = SkyCoord(ra=lon, dec=lat, frame='icrs', unit='deg')
        l, b = eq.galactic.l.value, eq.galactic.b.value
    else:
        l, b = lon, lat

    # conver to theta, phi
    theta = np.radians(90. - b)
    phi = np.radians(l)

    # convert to HEALPix indices
    indices = hp.ang2pix(nside, theta, phi)

    idx, counts = np.unique(indices, return_counts=True)

    # fill the fullsky map
    hpx_map = np.zeros(npix, dtype=int)
    hpx_map[idx] = counts

    return hpx_map

hpd_grid(sample, alpha=0.05, roundto=4)

Calculate highest posterior density (HPD) of array for given alpha. The HPD is the minimum width Bayesian credible interval (BCI). The function works for multimodal distributions, returning more than one mode

Parameters:

Name	Type	Description	Default
sample	Numpy array or python list	An array containing MCMC samples	required
alpha	float	Desired probability of type I error (defaults to 0.05)	0.05
roundto		Number of digits after the decimal point for the results	4

Returns:

Name	Type	Description
hpd	list with the highest density interval
x	array with grid points where the density was evaluated
y	array with the density values
modes	list listing the values of the modes

Source code in photXclus/utils.py

def hpd_grid(sample, alpha=0.05, roundto=4):
    """Calculate highest posterior density (HPD) of array for given alpha. 
    The HPD is the minimum width Bayesian credible interval (BCI). 
    The function works for multimodal distributions, returning more than one mode

    Parameters
    ----------

    sample : Numpy array or python list
        An array containing MCMC samples
    alpha : float
        Desired probability of type I error (defaults to 0.05)
    roundto: integer
        Number of digits after the decimal point for the results

    Returns
    ----------
    hpd: list with the highest density interval
    x: array with grid points where the density was evaluated
    y: array with the density values
    modes: list listing the values of the modes

    """
    sample = np.asarray(sample)
    sample = sample[~np.isnan(sample)]
    # get upper and lower bounds
    l = np.min(sample)
    u = np.max(sample)
    density = kde.gaussian_kde(sample)
    x = np.linspace(l, u, 2000)
    y = density.evaluate(x)
    #y = density.evaluate(x, l, u) waitting for PR to be accepted
    xy_zipped = zip(x, y/np.sum(y))
    xy = sorted(xy_zipped, key=lambda x: x[1], reverse=True)
    xy_cum_sum = 0
    hdv = []
    for val in xy:
        xy_cum_sum += val[1]
        hdv.append(val[0])
        if xy_cum_sum >= (1-alpha):
            break
    hdv.sort()
    #diff = (u-l)/20  # differences of 5%
    diff = (u-l)/100  # differences of 1%
    hpd = []
    hpd.append(round(min(hdv), roundto))
    for i in range(1, len(hdv)):
        if hdv[i]-hdv[i-1] >= diff:
            hpd.append(round(hdv[i-1], roundto))
            hpd.append(round(hdv[i], roundto))
    hpd.append(round(max(hdv), roundto))
    ite = iter(hpd)
    hpd = list(zip(ite, ite))
    modes = []
    for value in hpd:
         x_hpd = x[(x > value[0]) & (x < value[1])]
         y_hpd = y[(x > value[0]) & (x < value[1])]
         modes.append(round(x_hpd[np.argmax(y_hpd)], roundto))
    return hpd, x, y, modes