Commit bad06c61 authored by Akke Viitanen's avatar Akke Viitanen
Browse files

refactor software to src/

parent 9a8ef211

src/python/__init__.pye

deleted100644 → 0
+0 −0

Empty file deleted.

src/python/aird2018.pye

deleted100644 → 0
+0 −98
Original line number Diff line number Diff line 
#!/usr/bin/env python3

# encoding: utf-8

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2023-05-08 18:15:21



"""

Aird+2018 functions & data etc.

"""





import argparse

import glob

import math

import os

import random

import re

import subprocess

import sys

import time



import astropy as ap

import astropy.coordinates as c

import astropy.units as u

import fitsio

import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np

import scipy as sp



from util import ROOT





AIRD2018_TO = np.loadtxt(f"{ROOT}/data/aird2018/pledd_all.dat")

AIRD2018_QU = np.loadtxt(f"{ROOT}/data/aird2018/pledd_qu.dat")

AIRD2018_SF = np.loadtxt(f"{ROOT}/data/aird2018/pledd_sf.dat")





def get_aird2018(type):

    return {

        "all": AIRD2018_TO,

        "quiescent": AIRD2018_QU,

        "star-forming": AIRD2018_SF,

    }.get(type)





def get_zbins(type="all"):

    d = get_aird2018(type)

    zlo = np.unique(d[:, 0])

    zhi = np.unique(d[:, 1])

    return [(z1, z2) for z1, z2 in zip(zlo, zhi)]





def get_mbins(type="all"):

    d = get_aird2018(type)

    mlo = np.unique(d[:, 2])

    mhi = np.unique(d[:, 3])

    return [(m1, m2) for m1, m2 in zip(mlo, mhi)]





def get_l(type="all"):

    d = get_aird2018(type)

    return np.unique(d[:, 4])





def get_plambda(m, z, type="all"):



    air18 = get_aird2018(type)



    # NOTE: extrapolation to the nearest bin

    z = np.clip(z, air18[:, 0].min(), air18[:, 1].max() - 1e-9)

    select_z = (air18[:, 0] <= z) * (z < air18[:, 1])



    mmin = air18[select_z, 2].min() ; mmax = air18[select_z, 3].max()

    if   m <= mmin: m = mmin

    elif m >= mmax: m = mmax - 1e-9

    select_m = (air18[:, 2] <= m) * (m < air18[:, 3])



    # NOTE: returns [lambda, p(lambda)]

    return air18[select_z * select_m, 4:6].T





def get_duty_cycle(m, z, type="all"):

    dl = np.diff(get_l(type))

    assert np.allclose(dl[0], dl[1:])

    dl = dl[0]

    return np.sum(get_plambda(m, z, type) * dl)





def get_log_lambda_sBHAR(m, z, t):

    l_all = []

    for _m, _z, _t in zip(m, z, t):

        l_vec = get_l(_t)

        p = get_plambda(_m, _z, _t)

        c = np.cumsum(p) * np.diff(l_vec)[0]

        l = np.interp(np.random.rand(), c, l_vec, left=-np.inf, right=-np.inf)

        l_all.append(l)

    return np.array(l_all)

src/python/ananna2022.pye

deleted100644 → 0
+0 −215
Original line number Diff line number Diff line 
#!/usr/bin/env python3

# encoding: utf-8

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2024-08-01 13:53:22



"""

Implements Ananna+2022

"""





import argparse

import glob

import math

import os

import random

import re

import subprocess

import sys

import time



import astropy as ap

import astropy.coordinates as c

import astropy.units as u

import fitsio

import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np

import scipy as sp



labels = [

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

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

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

    r"Intrinsic ($\sigma=0.3; {\rm OA} = 35^{\circ}$)",

    r"$1/V_{\rm max}$",

]



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

# Ananna+ 2022, Table 3

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)

        ],



        "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)

        ],



        "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)

        ]



    },



    "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),



        ],



        "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),

        ],



        "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),

		]

    }



}



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



    """

    Return the Schechter function form of BHMF

    """



    if sample:

        

        ## NOTE: errors for sigma=0.50 case

        #m_star = 10 ** (np.log10(m_star) + np.random.uniform(-0.20, 0.25))

        #alpha += np.random.uniform(-0.110, +0.190)

        #beta  += np.random.uniform(-0.086, +0.065)



        m_star = 10 ** (np.log10(m_star) + np.random.normal(scale=0.50 * (0.20 + 0.25)))

        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

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



    # NOTE: fix for h

    #   original unit is 1/(Mpc/h)^3

    #   so that e.g. h = 0.70 corresponds to 1/(Mpc/0.70)^3 = 1/Mpc3 * 0.70^3



    return ret * h ** 3





def get_phi_lambda(l, l_star, xi_star, delta1, epsilon_lambda, h=1.0):

    ratio = l / l_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):



    x = np.log10(m)

    y = np.zeros_like(m)



    x1, y1 = np.loadtxt("data/ananna2022/fig10/bhmf/bhmf_ledd_gt_%d_type1.dat" % log_ledd_gt).T

    x2, y2 = np.loadtxt("data/ananna2022/fig10/bhmf/bhmf_ledd_gt_%d_type2.dat" % log_ledd_gt).T



    if is_type1: y += 10 ** np.interp(x, x1, y1, left=-np.inf, right=-np.inf)

    if is_type2: y += 10 ** np.interp(x, x2, y2, left=-np.inf, right=-np.inf)



    return y * h ** 3





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



    x = np.log10(l)

    y = np.zeros_like(l)



    x1, y1 = np.loadtxt("data/ananna2022/fig10/erdf/erdf_mbh_gt_%.1f_type1.dat" % log_mbh_gt).T

    x2, y2 = np.loadtxt("data/ananna2022/fig10/erdf/erdf_mbh_gt_%.1f_type2.dat" % log_mbh_gt).T



    if is_type1: y += 10 ** np.interp(x, x1, y1, left=-np.inf, right=-np.inf)

    if is_type2: y += 10 ** np.interp(x, x2, y2, left=-np.inf, right=-np.inf)



    return y * h ** 3





if __name__ == "__main__":



    m = np.logspace( 6.5,  10, 41)

    l = np.logspace(-3.0, 0.5, 36)

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



    for i in range(5):



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



            try:



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

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

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

                axes[0, j].set_xlim(10 ** 6.5, 10 ** 9.5)

                axes[0, j].set_ylim(1e-10, 1e-2)

                axes[0, j].set_title(k)

                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]

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

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

                axes[1, j].set_xlim(10 ** -3.0, 10 ** +0.5)

                axes[1, j].set_ylim(10 ** -8.0, 10 ** -2.5)

                axes[1, j].set_xlabel(r"$\lambda{\rm Edd}$")

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



            except IndexError as e:

                pass



    for ax in axes.flatten():

        ax.legend()



    plt.show()

    quit()



    plt.figure()

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

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

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

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

    plt.legend()

    plt.loglog()



    plt.figure()

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

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

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

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

    plt.legend()

    plt.loglog()

    plt.show()

src/python/baseline.py

+0 −1
Original line number Diff line number Diff line
@@ -43,4 +43,3 @@ def get_kwargs_instance_catalog(observation_id): 
        "numexp": 0,

        "seqnum": self.get_seqnum(observation_id),

    }

src/python/baseline.pye

deleted100644 → 0
+0 −46
Original line number Diff line number Diff line 
from astropy.table import Table

import sqlite3

import pandas as pd





def get_baseline(filename, query="SELECT * FROM observations", is_astropy=True):

    with sqlite3.open(f"{filename}") as con:

        ret = pd.read_sql(f"{query}", con)



    if is_astropy:

        ret = Table.from_pandas(red)



    return ret





def get_seqnum(observation_id, seqnum_max=32768):

    return observation_id % seqnum_max





def get_kwargs_instance_catalog(observation_id):

    kwargs_instance_catalog = {

        "mjd": d["observationStartMJD"],

        "rightascension": d["fieldRA"],

        "declination": d["fieldDec"],

        "altitude": d["altitude"],

        "azimuth": d["azimuth"],

        "filter": idx_filter,

        "rotskypos": d["rotSkyPos"],

        "dist2moon": d["moonDistance"],

        "moonalt": d["moonAlt"],

        "moondec": d["moonDec"],

        "moonphase": d["moonPhase"],

        "moonra": d["moonRA"],

        "nsnap": 1,

        #"obshistid": observationId,

        "obshistid": 0,

        "rottelpos": d["rotTelPos"],

        "seed": observation_id,

        #"seeing": d["seeingFwhm500"],

        "seeing": d["seeingFwhmEff"],

        "sunalt": d["sunAlt"],

        "vistime": d["visitExposureTime"],

        "numexp": 0,

        "seqnum": self.get_seqnum(observation_id),

    }