add schulze2015 regression tests

import numpy as np

from scipy.integrate import quad

def get_log_Psi(

    log_M_bh,

    log_lambda_edd,

    redshift,

    log_Psi_star=-5.32,

    log_M_bh_star=9.09,

    alpha=-1.50,

    beta=0.96,

    log_lambda_edd_star0=-1.19,

    alpha_lambda_edd=-0.29,

    k_lambda_edd=0.099,

    gamma=0.05,

    c_bh=0.270,

    c_lambda_edd=0.100,

    c_alpha_lambda_edd=0.094,

    log_M_bh_c=8,

    redshift_c=1.6,

):

    # Calculate log_rho_bh

    log_M_bh_star = log_M_bh_star + c_bh * (redshift - redshift_c)

    log_x = log_M_bh - log_M_bh_star

    log_rho_bh = (alpha + 1) * log_x - 10 ** (beta * log_x) / np.log(10) + np.log10(np.log(10))

    # Calculate log_rho_lambda_edd

    log_lambda_edd_star = (

        log_lambda_edd_star0 + k_lambda_edd * (log_M_bh - log_M_bh_c) + c_lambda_edd * (redshift - redshift_c)

    )

    # NOTE: typo in schulze? i.e. log_alpha_l should be alpha_l

    # log_alpha_l = alpha_l + c_alpha_l * (z - z_c)

    alpha_lambda_edd = alpha_lambda_edd + c_alpha_lambda_edd * (redshift - redshift_c)

    log_x = log_lambda_edd - log_lambda_edd_star

    log_rho_lambda_edd = (alpha_lambda_edd + 1) * log_x - 10**log_x / np.log(10) + np.log10(np.log(10))

    # Calculate log_rho_z

    log_rho_redshift = gamma * (redshift - redshift_c)

    return log_Psi_star + log_rho_bh + log_rho_lambda_edd + log_rho_redshift

@np.vectorize

def get_log_Phi_bh(log_M_bh, redshift, log_lambda_edd_min=-2, log_lambda_edd_max=1):

    def fun(log_lambda_edd):

        return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift)

    return np.log10(quad(fun, log_lambda_edd_min, log_lambda_edd_max)[0])

@np.vectorize

def get_log_Phi_lambda_edd(log_lambda_edd, redshift, log_M_bh_min=7, log_M_bh_max=11):

    def fun(log_M_bh):

        return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift)

    return np.log10(quad(fun, log_M_bh_min, log_M_bh_max)[0])

if __name__ == "__main__":

    log_m_bh, log_lambda_edd = np.meshgrid(np.linspace(7, 11), np.linspace(-2, 1))

    redshifts = 1.3, 1.9

    import matplotlib as mpl

    import matplotlib.pyplot as plt

    fig, axes = plt.subplots(1, 2, figsize=(2 * 6.4, 4.8))

    cmap = mpl.cm.get_cmap("viridis")

    norm = mpl.colors.BoundaryNorm([-11, -10, -9, -8, -7, -6, -5, -4], cmap.N)

    for i, redshift in enumerate(redshifts):

        log_Psi = get_log_Psi(log_m_bh, log_lambda_edd, redshift)

        im = axes[i].imshow(log_Psi, cmap=cmap, norm=norm, extent=(7, 11, -2, 1))

        plt.colorbar(im, ax=axes[i])

        axes[i].set_title("z = %.2f" % redshift)

    plt.savefig("fig/schulze2015_fig10.pdf")

import numpy as np

from lsst_inaf_agile import schulze2015

def test_get_log_Psi_one_value(num_regression):

    val = schulze2015.get_log_Psi(8.5, -0.5, 1.0)

    rec = {"result": val}

    num_regression.check(rec)

def test_get_log_Psi_2d(num_regression):

    x, y = np.meshgrid(np.linspace(8, 10), np.linspace(-2, 0))

    z = 1.2

    val = schulze2015.get_log_Psi(x, y, z).flatten()

    num_regression.check({"result": val})

def test_get_log_Psi_differet_redshifts(num_regression):

    val1 = schulze2015.get_log_Psi(8.5, -0.5, 1.0)

    val2 = schulze2015.get_log_Psi(8.5, -0.5, 1.2)

    val3 = schulze2015.get_log_Psi(8.5, -0.5, 1.4)

    num_regression.check(

        {

            "result1": val1,

            "result2": val2,

            "result3": val3,

        }

    )

def test_get_log_Psi_different_parameter_values(num_regression):

    val1 = schulze2015.get_log_Psi(8.5, -0.5, 1.0, -5.32)

    val2 = schulze2015.get_log_Psi(8.5, -0.5, 1.0, -5.22)

    val3 = schulze2015.get_log_Psi(8.5, -0.5, 1.0, -5.12)

    to_check = {

        "result1": val1,

        "result2": val2,

        "result3": val3,

    }

    np.random.seed(2022)

    for i in range(10):

        pvec = np.random.rand(13)

        val = schulze2015.get_log_Psi(8.5, -0.5, 1.0, *pvec)

        to_check[f"result_random{i}"] = val

    num_regression.check(to_check)

def test_get_log_Phi_bh(num_regression):

    args_all = [

        # change MBH

        (8.5, 1.0),

        (9.5, 1.0),

        (9.8, 1.0),

        # change lambda

        (8.5, 1.0),

        (8.5, 1.5),

        (8.5, 2.0),

        # change integration limits

        (8.5, 1.0, -2.0, 1.0),

        (8.5, 1.0, -2.0, 1.1),

        (8.5, 1.0, -2.1, 1.0),

        (8.5, 1.5, -2.0, 1.0),

        (9.0, 1.0, -2.0, 1.0),

    ]

    to_check = {

        f"result{i}": np.atleast_1d(schulze2015.get_log_Phi_bh(*args)) for i, args in enumerate(args_all)

    }

    num_regression.check(to_check)

def test_get_log_Phi_lambda_edd(num_regression):

    args_all = [

        # change MBH

        (-0.5, 1.0),

        (-0.3, 1.0),

        (-0.1, 1.0),

        # change lambda

        (-0.5, 1.0),

        (-0.5, 1.5),

        (-0.5, 2.0),

        # change integration limits

        (-0.5, 1.0, -2.0, 1.0),

        (-0.5, 1.0, -2.0, 1.1),

        (-0.5, 1.0, -2.1, 1.0),

        (-0.5, 1.5, -2.0, 1.0),

        (-0.5, 1.0, -2.0, 1.0),

    ]

    to_check = {

        f"result{i}": np.atleast_1d(schulze2015.get_log_Phi_lambda_edd(*args))

        for i, args in enumerate(args_all)

    }

    num_regression.check(to_check)

,result0,result1,result2,result3,result4,result5,result6,result7,result8,result9,result10

0,-4.9664629929571351,-6.8110538757359382,-8.4013722094776959,-4.9664629929571351,-4.8265002670457369,-4.6997983502920953,-4.9664629929571351,-4.9664629830534412,-4.9416287749002539,-4.8265002670457369,-5.5442273877711221

,result0,result1,result2,result3,result4,result5,result6,result7,result8,result9,result10

0,-5.9115519642393402,-7.2061785239129144,-9.1785460603467381,-5.9115519642393402,-5.5164700126971384,-5.1611066252051918,-12.486420069785996,-12.259674061910019,-12.486420069775743,-10.99045414348676,-12.486420069785996