implements zou f_occ and finish refactoring (dbc0a150) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

Makefile

+16 −0

Original line number	Diff line number	Diff line
		@@ -159,6 +159,9 @@ fig/flux_sigma_eta_ratio_20250707.pdf: src/python/plot_flux_sigma_eta_ratio.py
		fig/flux_sigma_eta_ratio_20250725.pdf: src/python/plot_flux_sigma_eta_ratio.py
		python3 $< $@ 0 \| tee table/table_flux_sigma_eta_ratio_20250725.tex

		fig/flux_sigma_eta_ratio_20250911.pdf: src/python/plot_flux_sigma_eta_ratio.py
		python3 $< $@ 0 \| tee table/table_flux_sigma_eta_ratio_20250911.tex

		fig/flux_sigma_eta_ratio_u.pdf: python/plot_flux_sigma_eta_ratio.py
		python3 $< $@ 0 u > table/flux_sigma_eta_u.table
		fig/flux_sigma_eta_ratio_g.pdf: python/plot_flux_sigma_eta_ratio.py
		@@ -238,3 +241,16 @@ diff:
		paper/version_sent_to_coauthors/manuscript_viitanen.tex \
		paper/version_after_coauthors/manuscript_viitanen.tex \
		> paper/version_after_coauthors/diff.tex

		test:
		PYTHONPATH=opt:src/python:$(PYTHONPATH) \
		python3 src/python/main.py --config etc/config_test.ini

		tests:
		PYTHONPATH=opt:src/python:$(PYTHONPATH)
		python3 src/tests/test_truth_catalog.py


		debug:
		PYTHONPATH=opt:src/python:$(PYTHONPATH) \
		python3 -m pdb src/python/main.py --config etc/config_test.ini

install.sh

+28 −2

Original line number	Diff line number	Diff line
		#!/bin/sh -ex
		#
		# Installation scripts for AGILE requirements
		#
		# Author: Akke Viitanen
		# Email: akke.viitanen@helsinki.fi

		@@ -38,6 +40,7 @@ install_vif () {
		fi
		}

		# Install EGG
		install_egg () {
		git clone git@github.com:cschreib/egg.git opt/egg
		cd opt/egg
		@@ -51,6 +54,7 @@ install_egg () {
		sudo chmod -Rv a+x /usr/local/share/egg
		}

		# Install imSim
		install_imsim () {

		# NOTE: only download the data
		@@ -62,14 +66,30 @@ install_imsim () {
		git clone https://github.com/LSSTDESC/imSim.git
		git clone https://github.com/LSSTDESC/skyCatalogs

		## conda dependencies:
		# mamba install -y --file imSim/etc/standalone_conda_requirements.txt
		# mamba install rubin-sim
		#
		## Install imSim:
		# pip install --no-deps imSim/
		# pip install --no-deps skyCatalogs/

		# install sims_sed_library
		mkdir -p rubin_sim_data/sims_sed_library
		curl 'https://s3df.slac.stanford.edu/groups/rubin/static/sim-data/rubin_sim_data/skybrightness_may_2021.tgz' \| tar -C rubin_sim_data -xz
		curl .https://s3df.slac.stanford.edu/groups/rubin/static/sim-data/rubin_sim_data/throughputs_2023_09_07.tgz. \| tar -C rubin_sim_data -xz
		curl .https://s3df.slac.stanford.edu/groups/rubin/static/sim-data/sed_library/seds_170124.tar.gz. \| tar -C rubin_sim_data/sims_sed_library -xz
		curl 'https://s3df.slac.stanford.edu/groups/rubin/static/sim-data/rubin_sim_data/throughputs_2023_09_07.tgz' \| tar -C rubin_sim_data -xz
		curl 'https://s3df.slac.stanford.edu/groups/rubin/static/sim-data/sed_library/seds_170124.tar.gz' \| tar -C rubin_sim_data/sims_sed_library -xz

		# setup envs
		conda env config vars set RUBIN_SIM_DATA_DIR=$(pwd)/rubin_sim_data
		conda env config vars set SIMS_SED_LIBRARY_DIR=$(pwd)/rubin_sim_data/sims_sed_library

		# finally, copy the const nu SED template
		cp ../../data/Constnu.Template.spec.gz rubin_sim_data/sims_sed_library

		}

		# Install LSST Science Pipelines
		install_lsst_science_pipelines () {

		# https://pipelines.lsst.io/install/lsstinstall.html
		@@ -86,6 +106,12 @@ install_lsst_science_pipelines () {
		eups distrib install -t v29_1_1 lsst_distrib
		curl -sSL 'https://raw.githubusercontent.com/lsst/shebangtron/main/shebangtron' \| python
		setup lsst_distrib

		# Download repo_assets
		# 1. create repo_assets
		# 2. download the files https://portal.nersc.gov/cfs/descssim/imSim/repo_assets/
		# 3. point REPO_ASSET_DIR to the location

		}

		#install_agile

src/python/catalog_agn.py

+74 −74

Original line number	Diff line number	Diff line
		@@ -51,7 +51,7 @@ class CatalogAGN:
		def __init__(
		self,
		dirname,
		egg,
		catalog_galaxy,
		type_plambda,
		save_sed,
		seed,
		@@ -60,9 +60,7 @@ class CatalogAGN:
		):

		self.dirname = dirname
		self.egg = egg
		self.bands = egg["BANDS"][0]
		self.rfbands = egg["RFBANDS"][0]
		self.catalog_galaxy = catalog_galaxy
		self.type_plambda = type_plambda
		self.save_sed = save_sed
		self.seed = seed
		@@ -79,7 +77,7 @@ class CatalogAGN:

		def __getitem__(self, key):
		if key not in self.catalog.dtype.names:
		return self.egg[key][0]
		return self.catalog_galaxy[key]
		return self.catalog[key]

		def _get_select(self, redshift=None, mstar=None, luminosity=None, t="all"):
		@@ -129,6 +127,7 @@ class CatalogAGN:
		("log_L_2500", np.float64, "Monochromatic 2500keV luminosity", "erg/s/Hz"),
		("MBH", np.float64, "Black hole mass", "Msun"),
		("log_L_bol", np.float64, "Bolometric AGN luminosity", "erg/s"),
		("occupation_fraction", np.float64, "SMBH occupation fraction", ""),
		] + [
		(b.strip() + "_point", np.float64, f"AGN {b} flux", "uJy") for b in bands
		] + [
		@@ -137,7 +136,10 @@ class CatalogAGN:

		def get_dtype(self):
		i = []
		for n, t, _, _ in self.get_columns(self.bands, self.rfbands):
		for n, t, _, _ in self.get_columns(
		self.catalog_galaxy.bands,
		self.catalog_galaxy.rfbands
		):
		i += [(n, t)]
		return np.dtype(i)

		@@ -148,10 +150,13 @@ class CatalogAGN:
		if os.path.exists(filename):
		return fitsio.read(filename)

		self.catalog = np.empty_like(self.egg["ID"][0], dtype=self.get_dtype())
		for col, _, _, _ in self.get_columns(self.bands, self.rfbands):
		if col in self.egg.columns.names:
		self.catalog[col] = self.egg[col][0]
		self.catalog = np.empty_like(self.catalog_galaxy["ID"], dtype=self.get_dtype())
		for col, _, _, _ in self.get_columns(
		self.catalog_galaxy.bands,
		self.catalog_galaxy.rfbands
		):
		if col in self.catalog_galaxy.get_dtype().names:
		self.catalog[col] = self.catalog_galaxy[col]
		else:
		self.catalog[col] = self.get_agn(col)

		@@ -171,7 +176,6 @@ class CatalogAGN:
		def get_agn(self, key):

		logger.info(f"Getting AGN attribute {key}")
		print(f"Getting AGN attribute {key}")

		ret = None
		func = getattr(self, f"get_{key}", None)
		@@ -201,38 +205,37 @@ class CatalogAGN:
		if idxs is None:
		idxs = np.arange(self["ID"].size)

		log_lambda_SAR = np.full_like(self["Z"], np.nan)
		# Get the relevant galaxy properties
		m = self["M"]
		z = self["Z"]
		p = self["PASSIVE"]
		log_lambda_SAR = np.full_like(z, np.nan)

		# NOTE: short-circuit for CTN only as it is significantly faster
		if not add_ctk:
		logger.info("Assigning log_lambda_SAR for CTN...")
		for t in ["star-forming", "quiescent"]:
		select = self._get_select(t=t)
		U = np.random.rand(select.sum())
		log_lambda_SAR[select] = zou2024.get_inv_log_Plambda(
		np.log10(U),
		m[select],
		z[select],
		t
		np.log10(U), m[select], z[select], t
		)
		else:
		return log_lambda_SAR

		logger.info("Building the arguments...")
		arguments = [
		(i, self.catalog.size, self["M"][i], self["Z"][i], "quiescent" if self["PASSIVE"][i] == 1 else "star-forming")
		(i, self.catalog.size, m[i], z[i], "quiescent" if p[i] == 1 else "star-forming")
		for i in np.arange(self.catalog.size)
		]

		cpu_count = multiprocessing.cpu_count()
		logger.info(f"Assigning log_lambda_SAR with {cpu_count=}...")
		logger.info(f"Assigning log_lambda_SAR N={log_lambda_SAR.size} with {cpu_count=}...")

		t0 = time.time()
		with multiprocessing.Pool(processes=cpu_count) as p:
		log_lambda_SAR = p.starmap(util.get_log_lambda_SAR, arguments)

		t1 = time.time()

		logger.info(f"assigned {len(log_lambda_SAR)} objects in {t1 - t0} seconds")
		logger.info(f"Assigned {len(log_lambda_SAR)} objects in {t1 - t0} seconds")

		return log_lambda_SAR

		@@ -246,9 +249,10 @@ class CatalogAGN:
		return x + self["M"]

		def get_log_FX_2_10(self, Gamma=1.9):
		# NOTE: 1 + z dependence from K-correction
		z = self["Z"].astype(np.float64)
		d = self["D"].astype(np.float64) * u.Mpc.to(u.cm)

		# NOTE: 1 + z dependence from K-correction
		return self["log_LX_2_10"] - np.log10(4 * np.pi * d ** 2) + (2 - Gamma) * np.log10(1 + z)

		def get_log_FX_2_7(self, Gamma=1.9):
		@@ -434,11 +438,6 @@ class CatalogAGN:
		wavlen = np.append(wavlen, [1450, 4400, 5007, 150000])
		wavlen = np.sort(wavlen)

		# NOTE: we need the egg 1000-4000 flux for the host galaxy test. We
		# recover the column name here for efficiency
		idx_mock_1000_4000 = np.where(self.egg["RFBANDS"][0] == "mock-1000-4000")[0][0]
		assert idx_mock_1000_4000 >= 0

		# Populate the SED
		while True:

		@@ -454,7 +453,7 @@ class CatalogAGN:
		LogL2kev=self["log_L_2_keV"][i],
		add_NL=self["is_optical_type2"][i],
		NL_normalization="lamastra",
		Av_lines=self.egg["AVLINES_BULGE"][0, i] + self.egg["AVLINES_DISK"][0, i],
		Av_lines=self.catalog_galaxy["AVLINES_BULGE"][i] + self.catalog_galaxy["AVLINES_DISK"][i],
		*dict(zip(PARAMETER_NAMES, POSTERIOR_DISTRIBUTION.sample().values))
		)

		@@ -473,7 +472,7 @@ class CatalogAGN:
		band="mock-1000-4000",
		filter_db=self.filter_db
		)
		mag_gal = self.egg["RFMAG"][0, i, idx_mock_1000_4000]
		mag_gal = self.catalog_galaxy["magabs_mock-1000-4000"][i]
		flux_gal = util.mag_to_flux(mag_gal)
		ratio = flux_agn / (flux_gal + flux_agn)

		@@ -504,13 +503,9 @@ class CatalogAGN:
		sed.write(filename, lam, flux)

		# NOTE: optimization: save values from the sed for future use
		for b, r in product(map(str.strip, self.bands), [False, True]):
		for b, r in product(map(str.strip, self.catalog_galaxy.bands), [False, True]):
		FLUXES[i, b, r] = self._get_flux_single(
		agn_sed,
		b,
		r,
		self["Z"][i],
		self["D"][i]
		agn_sed, b, r, self["Z"][i], self["D"][i]
		)

		return agn_sed
		@@ -575,22 +570,11 @@ class CatalogAGN:
		select=None
		):

		# Get the selection
		if select is None:
		#select = np.ones_like(self["M"], dtype=bool)
		select = np.isfinite(self["Z"])

		# Get MBH
		log_mbh = np.full_like(self["M"], np.nan)

		# 20250508 updated p(lambda) Mbh-Mstar
		#log_mbh[select] = np.array([get_log_mbh_continuity(c["M"], c["Z"]) for c in self.catalog[select]])
		#log_mbh[select] = get_log_mbh_continuity_new(self["M"][select], self["Z"][select])
		log_mbh[select] = get_log_mbh_continuity_new2(self["M"][select], self["Z"][select])
		log_mbh = get_log_mbh_continuity_new2(self["M"], self["Z"])

		# Add the scatter
		sigma_tot = np.full_like(self["M"], -np.inf)
		sigma_tot[select] = np.random.normal(scale=sigma_total, size=select.sum())
		sigma_tot = np.random.normal(scale=sigma_total, size=log_mbh.size)

		# NOTE: add the offset?
		return log_mbh + sigma_tot + offset
		@@ -643,6 +627,23 @@ class CatalogAGN:
		a = 12.6236
		return self[key] - a * self["E_BV"]

		def get_occupation_fraction(self):
		"""
		Returns the black hole occupation fraction 'f_occ' as a function of
		logMstar. The occupation fraction is the percentage of host galaxies
		that are expected to host SMBH, which could be <1 for dwarf galaxies
		(Miller+15, Burke+25, Zou+).

		The updated AGN duty cycle can be estimated using:

		is_agn_new = occupation_fraction(logMstar) * is_agn_old
		"""

		table = Table.read("data/from_zou/table2.csv")
		x = table["logm"]
		y = table["median"]
		return np.interp(self["M"], x, y, left=0.0, right=1.0)

		def get_lightcurve(
		self,
		i,
		@@ -678,7 +679,6 @@ class CatalogAGN:
		np.save(filename, lc)
		return lc


		def get_lightcurve_mjd(self, i, band, mjd, mjd0):
		lc = self.get_lightcurve(i, band)
		return np.interp(mjd - mjd0, util.get_mjd_vec(), lc)

src/python/catalog_combined.py

+6 −8

Original line number	Diff line number	Diff line
		@@ -70,7 +70,6 @@ class CatalogCombined:
		len(self.catalog_binary.catalog)
		)


		def get_catalog_combined(self):

		"""
		@@ -191,6 +190,5 @@ class CatalogCombined:
		# Handles something else
		return None


		def __getitem__(self, key):
		return self.catalog_combined[key]

src/python/catalog_galaxy.py

+44 −11

Original line number	Diff line number	Diff line
		@@ -19,7 +19,15 @@ logger = logging.getLogger(__name__)

		class CatalogGalaxy:

		def __init__(self, dirname, egg):
		"""
		Galaxy catalog
		"""

		def __init__(self, dirname: str, egg: dict):

		"""
		Initialize the galaxy catalog
		"""

		self.dirname = dirname
		self.egg = egg
		@@ -30,6 +38,10 @@ class CatalogGalaxy:

		def get_catalog(self):

		"""
		Generates the galaxy catalog and writes it to disk
		"""

		# Create the Galaxy catalog
		self.catalog = np.empty_like(self.egg["RA"][0], dtype=self.get_dtype())
		filename = f"{self.dirname}/galaxy.fits"
		@@ -51,12 +63,18 @@ class CatalogGalaxy:
		@staticmethod
		def get_columns(bands, rfbands, is_public=False):

		"""
		Returns the galaxy catalog columns, types, and descriptions. The
		arguments 'bands' and 'rfbands' are used to generate the columns for
		the fluxes and apparent magnitudes, respectively.
		"""

		ret = [
		("ID", np.int64, "unique ID", ""),
		("RA", np.float64, "Right ascenscion", "deg"),
		("DEC", np.float64, "Declination", "deg"),
		("Z", np.float64, "Cosmological redshift", ""),
		("D", np.float64, "Luminosity distance OR distance modulus for stars", "Mpc OR mag"),
		("D", np.float64, "Luminosity distance", "Mpc"),
		("M", np.float64, "log10 of stellar mass", "Msun"),
		("SFR", np.float64, "Star-formation rate", "Msun/yr"),
		("PASSIVE", bool, "Is passive (non-star-forming)", ""),
		@@ -79,6 +97,11 @@ class CatalogGalaxy:
		return ret

		def get_dtype(self):

		"""
		Returns the numpy dtype corresponding to the columns
		"""

		dtype = []
		for n, t, _, _ in self.get_columns(self.bands, self.rfbands):
		dtype += [(n, t)]
		@@ -86,6 +109,10 @@ class CatalogGalaxy:

		def get_galaxy(self, key):

		"""
		Returns the galaxy property corresponding to 'key'
		"""

		logger.info(f"Getting galaxy attribute {key}")

		_bands = [b.strip() for b in self.bands]
		@@ -112,11 +139,12 @@ class CatalogGalaxy:
		# Fluxes not generated... return zerovalues
		return np.full_like(self.egg["Z"][0], np.inf if "magabs_" in key else 0.0)

		ret = self.egg[key][0]
		# Handles stellar mass
		if key == "M":
		# NOTE: Salpeter to Chabrier IMF conversion see e.g. Grylls+20
		ret -= 0.24
		return self.egg["M"][0] - 0.24

		# Handles star-formation rate
		if key == "SFR":
		# Where necessary to convert SFRs from the literature from Chabrier
		# or Kroupa IMFs to the Salpeter IMF, we divide by constant factors
		@@ -125,8 +153,13 @@ class CatalogGalaxy:
		# SFR_salpeter = SFR_chabrier / 0.63
		# ==> SFR_chabrier = 0.63 * SFR_salpeter
		# (https://ned.ipac.caltech.edu/level5/March14/Madau/Madau3.html)
		ret *= 0.63
		return ret
		return self.egg["SFR"][0] * 0.63

		# Handles any other egg column
		return self.egg[key][0]

		def __getitem__(self, key):
		"""
		Returns the galaxy catalog column corresponding to 'key'
		"""
		return self.get_galaxy(key)