Unverified Commit c89ad561 authored by Akke Viitanen's avatar Akke Viitanen
Browse files

adds unit tests

parent fdc2bcfc

Makefile

+3 −0
Original line number Diff line number Diff line
@@ -270,3 +270,6 @@ figures_20250930: 
	python3 src/scripts/plots/plot_black_hole_mass_function.py

	# gband number counts

	python3 src/scripts/plots/plot_number_counts_gband.py



coverage:

	pytest --cov=src --cov-report=html tests/

src/lsst_inaf_agile/ananna2022.py

+64 −40
Original line number Diff line number Diff line
@@ -7,8 +7,9 @@ 


import matplotlib.pyplot as plt

import numpy as np

from numpy.typing import ArrayLike



labels = [

LABELS = [

    r"Intrinsic ($\sigma=0.3$)",

    r"Intrinsic ($\sigma=0.3; \sigma_{\log L,{\rm scatt}} = 0.2$)",

    r"Intrinsic ($\sigma=0.5$)",
@@ -18,55 +19,63 @@ labels = [ 


###############################################################################

# Ananna+ 2022, Table 3

parameters = {

PARAMETERS = {

    "BHMF": {

        "All": [

            (labels[0], 10**7.88, 10**-3.52, -1.576, 0.593),

            (labels[1], 10**7.92, 10**-3.67, -1.530, 0.612),

            (labels[2], 10**7.67, 10**-3.37, -1.260, 0.630),

            (labels[3], 10**7.92, 10**-3.49, -1.576, 0.600),

            (labels[4], 10**8.12, 10**-4.33, -1.060, 0.574),

            (LABELS[0], 10**7.88, 10**-3.52, -1.576, 0.593),

            (LABELS[1], 10**7.92, 10**-3.67, -1.530, 0.612),

            (LABELS[2], 10**7.67, 10**-3.37, -1.260, 0.630),

            (LABELS[3], 10**7.92, 10**-3.49, -1.576, 0.600),

            (LABELS[4], 10**8.12, 10**-4.33, -1.060, 0.574),

        ],

        "Type 1": [

            (labels[0], 10**7.97, 10**-4.19, -1.753, 0.561),

            (labels[1], 10**7.93, 10**-4.27, -1.730, 0.566),

            (labels[2], 10**7.91, 10**-4.27, -1.560, 0.590),

            (labels[4], 10**8.73, 10**-5.10, -1.350, 0.681),

            (LABELS[0], 10**7.97, 10**-4.19, -1.753, 0.561),

            (LABELS[1], 10**7.93, 10**-4.27, -1.730, 0.566),

            (LABELS[2], 10**7.91, 10**-4.27, -1.560, 0.590),

            (LABELS[4], 10**8.73, 10**-5.10, -1.350, 0.681),

        ],

        "Type 2": [

            (labels[0], 10**7.820, 10**-3.60, -1.16, 0.637),

            (labels[1], 10**7.790, 10**-3.64, -1.18, 0.617),

            (labels[2], 10**7.760, 10**-3.60, -0.99, 0.703),

            (labels[3], 10**7.730, 10**-3.44, -1.26, 0.635),

            (labels[4], 10**8.102, 10**-4.33, -1.04, 0.732),

            (LABELS[0], 10**7.820, 10**-3.60, -1.16, 0.637),

            (LABELS[1], 10**7.790, 10**-3.64, -1.18, 0.617),

            (LABELS[2], 10**7.760, 10**-3.60, -0.99, 0.703),

            (LABELS[3], 10**7.730, 10**-3.44, -1.26, 0.635),

            (LABELS[4], 10**8.102, 10**-4.33, -1.04, 0.732),

        ],

    },

    "ERDF": {

        "All": [

            (labels[0], 10**-1.338, 10**-3.64, 0.38, 2.260),

            (labels[1], 10**-1.286, 10**-3.76, 0.40, 2.322),

            (labels[2], 10**-1.332, 10**-3.68, 0.484, 2.210),

            (labels[3], 10**-1.249, 10**-3.80, 0.28, 2.720),

            (labels[4], 10**-1.190, 10**-3.76, -0.02, 2.060),

            (LABELS[0], 10**-1.338, 10**-3.64, 0.38, 2.260),

            (LABELS[1], 10**-1.286, 10**-3.76, 0.40, 2.322),

            (LABELS[2], 10**-1.332, 10**-3.68, 0.484, 2.210),

            (LABELS[3], 10**-1.249, 10**-3.80, 0.28, 2.720),

            (LABELS[4], 10**-1.190, 10**-3.76, -0.02, 2.060),

        ],

        "Type 1": [

            (labels[0], 10**-1.152, 10**-4.08, 0.30, 2.51),

            (labels[1], 10**-1.138, 10**-4.09, 0.27, 2.57),

            (labels[2], 10**-1.103, 10**-4.23, 0.13, 2.97),

            (labels[4], 10**-1.060, 10**-4.02, -0.51, 2.57),

            (LABELS[0], 10**-1.152, 10**-4.08, 0.30, 2.51),

            (LABELS[1], 10**-1.138, 10**-4.09, 0.27, 2.57),

            (LABELS[2], 10**-1.103, 10**-4.23, 0.13, 2.97),

            (LABELS[4], 10**-1.060, 10**-4.02, -0.51, 2.57),

        ],

        "Type 2": [

            (labels[0], 10**-1.657, 10**-3.82, 0.376, 2.50),

            (labels[1], 10**-1.628, 10**-3.84, 0.320, 2.50),

            (labels[2], 10**-1.675, 10**-3.80, 0.330, 2.51),

            (labels[3], 10**-1.593, 10**-3.92, 0.300, 2.53),

            (labels[4], 10**-1.870, 10**-3.74, -0.500, 2.30),

            (LABELS[0], 10**-1.657, 10**-3.82, 0.376, 2.50),

            (LABELS[1], 10**-1.628, 10**-3.84, 0.320, 2.50),

            (LABELS[2], 10**-1.675, 10**-3.80, 0.330, 2.51),

            (LABELS[3], 10**-1.593, 10**-3.92, 0.300, 2.53),

            (LABELS[4], 10**-1.870, 10**-3.74, -0.500, 2.30),

        ],

    },

}





def get_phi_bh(m, m_star, phi_star, alpha, beta, h=1.0, sample=None):

def get_phi_bh(

    m: ArrayLike,

    m_star: float,

    phi_star: float,

    alpha: float,

    beta: float,

    h: float = 1.0,

    sample: bool = False,

) -> ArrayLike:

    """Return the Schechter function form of BHMF."""

    if sample:

        ## NOTE: errors for sigma=0.50 case
@@ -78,7 +87,7 @@ def get_phi_bh(m, m_star, phi_star, alpha, beta, h=1.0, sample=None): 
        alpha += np.random.normal(scale=0.50 * (0.110 + 0.190))

        beta += np.random.normal(scale=0.50 * (0.086 + 0.065))



    x = m / m_star

    x = np.asarray(m) / m_star

    ret = np.log(10) * phi_star * x ** (alpha + 1) * np.exp(-(x**beta))



    # NOTE: fix for h
@@ -88,15 +97,24 @@ def get_phi_bh(m, m_star, phi_star, alpha, beta, h=1.0, sample=None): 
    return ret * h**3





def get_phi_lambda(lambda_edd, lambda_edd_star, xi_star, delta1, epsilon_lambda, h=1.0):

def get_phi_lambda(

    lambda_edd: ArrayLike,

    lambda_edd_star: float,

    xi_star: float,

    delta1: float,

    epsilon_lambda: float,

    h: float = 1.0,

) -> ArrayLike:

    """Return phi_lambda following the functional form in Ananna."""

    ratio = lambda_edd / lambda_edd_star

    ratio = np.asarray(lambda_edd) / lambda_edd_star

    return (

        np.ma.true_divide(xi_star, np.power(ratio, delta1) + np.power(ratio, delta1 + epsilon_lambda)) * h**3

    )





def get_phi_bh_fig10(m, is_type1=True, is_type2=True, log_ledd_gt=-3, h=1.0):

def get_phi_bh_fig10(

    m: ArrayLike, is_type1: bool = True, is_type2: bool = True, log_ledd_gt: float = -3, h: float = 1.0

):

    """Get phi_bh from Ananna Fig10."""

    x = np.log10(m)

    y = np.zeros_like(m)
@@ -112,7 +130,13 @@ def get_phi_bh_fig10(m, is_type1=True, is_type2=True, log_ledd_gt=-3, h=1.0): 
    return y * h**3





def get_phi_lambda_fig10(lambda_edd, is_type1=True, is_type2=True, log_mbh_gt=6.5, h=1.0):

def get_phi_lambda_fig10(

    lambda_edd: ArrayLike,

    is_type1: bool = True,

    is_type2: bool = True,

    log_mbh_gt: float = 6.5,

    h: float = 1.0,

) -> ArrayLike:

    """Get phi_lambda from Ananna Fig10."""

    x = np.log10(lambda_edd)

    y = np.zeros_like(lambda_edd)
@@ -136,7 +160,7 @@ if __name__ == "__main__": 
    for i in range(5):

        for j, k in enumerate(["All", "Type 1", "Type 2"]):

            try:

                p = parameters["BHMF"][k][i]

                p = PARAMETERS["BHMF"][k][i]

                phi = get_phi_bh(mbh, *p[1:])

                axes[0, j].loglog(mbh, phi, label=p[0])

                axes[0, j].set_xlim(10**6.5, 10**9.5)
@@ -145,7 +169,7 @@ if __name__ == "__main__": 
                axes[0, j].set_xlabel(r"$M_{\rm BH}$ [Msun]")

                axes[0, j].set_ylabel(r"$\Phi_{\rm BH}$ [1/(Mpc3/h3)/dex")



                p = parameters["ERDF"][k][i]

                p = PARAMETERS["ERDF"][k][i]

                phi = get_phi_lambda(lambda_edd, *p[1:])

                axes[1, j].loglog(lambda_edd, phi, label=p[0])

                axes[1, j].set_xlim(10**-3.0, 10**+0.5)
@@ -163,7 +187,7 @@ if __name__ == "__main__": 
    quit()



    plt.figure()

    for p in parameters["BHMF"]["All"]:

    for p in PARAMETERS["BHMF"]["All"]:

        plt.plot(mbh, get_phi_bh(mbh, *p[1:]), label=p[0])

    plt.xlabel(r"$M_{\rm BH}$ [Msun]")

    plt.ylabel(r"$\Phi_{M}$ [1/(Mpc/h)$^3$/dex]")
@@ -171,7 +195,7 @@ if __name__ == "__main__": 
    plt.loglog()



    plt.figure()

    for p in parameters["ERDF"]["All"]:

    for p in PARAMETERS["ERDF"]["All"]:

        plt.plot(lambda_edd, get_phi_lambda(lambda_edd, *p[1:]), label=p[0])

    plt.xlabel(r"$\lambda$")

    plt.ylabel(r"$\Phi_{\lambda}$ [1/(Mpc/h)$^3$/dex]")