Separates different catalog classes and image_simu (83f121dd) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

src/python/catalog_agn.py

+65 −95

Original line number	Diff line number	Diff line
		@@ -32,13 +32,12 @@ import zou2024

		logger = logging.getLogger(__name__)

		# Optimization: store the pre-computed AGN fluxes and lightcurves
		# Optimization: store the pre-computed AGN fluxes
		FLUXES = {}
		LIGHTCURVES = {}

		try:
		# Pre-load the posterior distribution for the AGN SED
		#logger.info("Reading in the posterior distribution...")
		logger.info("Reading in the posterior distribution...")
		#POSTERIOR_DISTRIBUTION = pd.read_csv(f"{ROOT}/posteriors/posterior.dat", sep=' ')
		#POSTERIOR_DISTRIBUTION = pd.read_csv(f"{ROOT}/data/posteriors/posterior_2024_05_20.dat", sep=' ')
		POSTERIOR_DISTRIBUTION = pd.read_csv(f"{ROOT}/data/posteriors/posterior_frozen.dat", sep=' ')
		@@ -49,30 +48,35 @@ except FileNotFoundError:

		class CatalogAGN:

		def __init__(self, dirname, egg, type_plambda, save_sed, seed, merloni2014):
		def __init__(
		self,
		dirname,
		egg,
		type_plambda,
		save_sed,
		seed,
		merloni2014,
		filter_db="data/egg/share/filter-db/db.dat"
		):

		self.dirname = dirname
		self.egg = egg
		self.bands = egg["BANDS"][0]
		self.rfbands = egg["RFBANDS"][0]
		self.type_plambda = type_plambda
		self.save_sed = int(save_sed)
		self.save_sed = save_sed
		self.seed = seed
		self.merloni2014 = merloni2014
		self.filter_db = filter_db

		# Set the seed
		if self.seed is not None:
		self.seed = int(self.seed)
		np.random.seed(self.seed)


		# Create the catalog
		self.catalog = self.get_catalog()

		# Create the lightcurves
		for b in "ugrizy":
		self.get_lightcurve_agn(band=f"lsst-{b}")

		def __getitem__(self, key):
		if key not in self.catalog.dtype.names:
		return self.egg[key][0]
		@@ -465,7 +469,12 @@ class CatalogAGN:
		redshift=0.0,
		distance_in_cm=10 * u.pc.to(u.cm),
		)
		flux_agn = sed.get_flux_band(lam, flux_agn, band="mock-1000-4000")
		flux_agn = sed.get_flux_band(
		lam,
		flux_agn,
		band="mock-1000-4000",
		filter_db=self.filter_db
		)
		mag_gal = self.egg["RFMAG"][0, i, idx_mock_1000_4000]
		flux_gal = util.mag_to_flux(mag_gal)
		ratio = flux_agn / (flux_gal + flux_agn)
		@@ -498,7 +507,7 @@ class CatalogAGN:

		# NOTE: optimization: save values from the sed for future use
		for b, r in product(map(str.strip, self.bands), [False, True]):
		FLUXES[i, b, r] = self._init_flux_single(
		FLUXES[i, b, r] = self._get_flux_single(
		agn_sed,
		b,
		r,
		@@ -508,7 +517,7 @@ class CatalogAGN:

		return agn_sed

		def _init_flux_single(self, my_sed, band, rest_frame, redshift, distance):
		def _get_flux_single(self, my_sed, band, rest_frame, redshift, distance):

		# Get the flux in observed frame or rest frame
		if rest_frame:
		@@ -523,7 +532,9 @@ class CatalogAGN:
		)

		# Get the flux for the requested band
		flux_band = sed.get_flux_band(lam, flux, band)
		flux_band = sed.get_flux_band(
		lam, flux, band, filter_db=self.filter_db
		)

		# NOTE: convert to magnitudes for the rest_frame
		if rest_frame:
		@@ -557,87 +568,6 @@ class CatalogAGN:

		return log_L_bol

		def _get_lightcurve_agn(self, i, band, args, *kwargs):

		kwargs.update({
		"mjd0": 0,
		"mjd": np.arange(0, 3653, 1),
		"band": band,
		"flux": self[f"{band}_point"][i],
		"z": self["Z"][i],
		# 20241004:
		#"mag_i": self["magabs_lsst-i"][i],
		"mag_i": self["magabs_lsst-i_point"][i],
		"logMbh": self["MBH"][i],
		"type2": self["is_optical_type2"][i],
		"T": 3653,
		"deltatc": 1,
		"seed": i + 1
		})

		logger.debug(f"Estimating AGN lightcurve {i} {band}")

		lc = lightcurve.get_lightcurve_agn(args, *kwargs)

		return i, band, lc

		def get_lightcurve_agn(
		self,
		idxs=None,
		band="lsst-r",
		T=3653,
		deltatc=1,
		mjd0=0,
		mjd=None,
		use_existing=True, # use an existing lightcurve
		overwrite=False, # overwrite existing lightcurve
		*args,
		**kwargs
		):

		"""
		Return the AGN lightcurves for the requested uids
		"""

		if mjd is None:
		mjd = np.arange(mjd0, mjd0 + T, deltatc)

		# NOTE: select all AGN if idxs is None
		if idxs is None:
		is_agn = self["is_agn"]
		uid_agn = self["ID"][is_agn]
		idxs = uid_agn
		idxs = np.atleast_1d(idxs)

		logger.debug(f"Getting AGN lightcurves {band} {idxs} {idxs.size}")

		for ndone, i in enumerate(idxs):

		logger.debug(f"{i} {ndone} {idxs.size} {100 * ndone / idxs.size:7.2}")
		filename = f"{self.dirname}/lightcurves/{i}/{band}.npy"

		if not os.path.exists(os.path.dirname(filename)):
		os.makedirs(os.path.dirname(filename))

		lc = None
		if not use_existing:
		_, _, lc = self._get_lightcurve_agn(i, band)
		if overwrite and os.path.exists(filename):
		np.save(filename, lc)
		# 20241029 patch...
		elif not os.path.exists(filename):
		_, _, lc = self._get_lightcurve_agn(i, band)
		np.save(filename, lc)
		else:
		lc = np.load(filename)

		LIGHTCURVES[i, band] = lc

		# Initialize the time and return the interpolated flux
		tt = np.arange(0, T, deltatc)
		ret = np.array([np.interp(mjd, tt + mjd0, LIGHTCURVES[i, band]) for i in idxs])
		return mjd, ret

		def get_MBH(
		self,
		sigma_total=0.50,
		@@ -714,3 +644,43 @@ class CatalogAGN:
		elif "mock-1450" in key:
		a = 12.6236
		return self[key] - a * self["E_BV"]

		def get_lightcurve(
		self,
		i,
		band,
		mjd=util.get_mjd_vec(),
		*args,
		**kwargs
		):

		filename = f"{self.dirname}/lightcurves/agn/{i}/{band}.npy"
		if os.path.exists(filename):
		logger.info(f"Reading AGN lightcurve {filename}")
		return np.load(filename)

		kwargs.update({
		"mjd0": 0.0,
		"mjd": mjd,
		"band": band,
		"flux": self[f"{band}_point"][i],
		"z": self["Z"][i],
		"mag_i": self["magabs_lsst-i_point"][i],
		"logMbh": self["MBH"][i],
		"type2": self["is_optical_type2"][i],
		"T": mjd.max() + 1,
		"deltatc": 1,
		"seed": self._get_seed("get_lightcurve") + i
		})
		logger.info(f"Estimating AGN lightcurve {i} {band}")
		lc = lightcurve.get_lightcurve_agn(args, *kwargs)

		logger.info(f"Writing AGN lightcurve {filename}")
		util.create_directory(filename)
		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_binary.py

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python3
		# encoding: utf-8
		# Author: Akke Viitanen
		# Email: akke.viitanen@helsinki.fi
		# Date: 2025-07-03 19:48:27

		"""
		"""

src/python/catalog_combined.py

0 → 100644

+189 −0

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-02-13 16:20:51

		"""
		Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution
		"""

		import logging

		import fitsio
		import numpy as np

		import util

		logger = logging.getLogger(__name__)


		class CatalogCombined:

		def __init__(self, dirname, catalog_agn, catalog_galaxy, catalog_star, catalog_binary):

		self.dirname = dirname
		self.filename = f"{dirname}/catalog.fits"

		self.catalog_agn = catalog_agn
		self.catalog_galaxy = catalog_galaxy
		self.catalog_star = catalog_star
		self.catalog_binary = catalog_binary

		self.catalog_combined = self.get_catalog_combined()
		self.catalog_combined = self.postprocess()
		self.write()

		def get_dtype(self):
		dtype = {}
		for catalog in (
		self.catalog_galaxy,
		self.catalog_agn,
		self.catalog_star,
		self.catalog_binary
		):
		for name in catalog.catalog.dtype.names:
		dtype[name] = catalog.catalog.dtype[name]

		for band in self.catalog_galaxy.bands:
		dtype[band + "_total"] = np.float64

		for band in self.catalog_galaxy.rfbands:
		dtype["magabs_" + band + "_total"] = np.float64

		return [(k, v) for k, v in dtype.items()]

		def get_number_galaxy_star_binary(self):
		"""
		Return the number of galaxies, stars, and binaries
		"""
		return (
		len(self.catalog_galaxy.catalog),
		len(self.catalog_star.catalog),
		len(self.catalog_binary.catalog)
		)


		def get_catalog_combined(self):

		"""
		Constructs the combined catalog based on the dtype.

		For each column in the catalog, the column is looked for in the
		AGN/galaxy/star/binary catalogs. Upon a match, the corresponding rows
		are assigned the values from the reference catalog. Note that some
		columns may appear in multiple catalog.

		The indices correspond to the following

		type idx_first idx_last

		AGN 0 n_galaxy
		galaxy 0 n_galaxy
		star n_galaxy n_galaxy + n_star
		binary n_galaxy + n_star n_total
		"""

		n_galaxy, n_star, n_binary = self.get_number_galaxy_star_binary()
		n_total = n_galaxy + n_star + n_binary
		catalog_combined = np.full(n_total, np.nan, dtype=self.get_dtype())

		# Set the catalog ID
		catalog_combined["ID"] = np.arange(n_total)

		# Fetch the values for the columns
		for name in catalog_combined.dtype.names:

		logger.info(f"Setting attribute {name} {catalog_combined.dtype[name]}")

		# Handles galaxy
		if name in self.catalog_agn.catalog.dtype.names:
		catalog_combined[name][:n_galaxy] = self.catalog_agn[name]

		# Handles AGN
		if name in self.catalog_galaxy.catalog.dtype.names:
		catalog_combined[name][:n_galaxy] = self.catalog_galaxy[name]

		# Handles star
		if name in self.catalog_star.catalog.dtype.names:
		catalog_combined[name][n_galaxy:n_galaxy + n_star] = self.catalog_star[name]
		elif catalog_combined.dtype[name] == np.dtype('bool'):
		catalog_combined[name][n_galaxy:n_galaxy + n_star] = False

		# Handles binary star
		if name in self.catalog_binary.catalog.dtype.names:
		catalog_combined[name][n_galaxy + n_star:] = self.catalog_binary[name]
		elif catalog_combined.dtype[name] == np.dtype('bool'):
		catalog_combined[name][n_galaxy + n_star:] = False

		return catalog_combined

		def postprocess(self):

		"""
		Adds columns in post-processing
		"""

		for band in self.catalog_galaxy.bands:
		logger.info(f"Postprocessing flux {band}")
		self.catalog_combined[band + "_total"] = self.get_flux_total(band, False)

		for band in self.catalog_galaxy.rfbands:
		logger.info(f"Postprocessing magabs {band}")
		self.catalog_combined["magabs_" + band + "_total"] = self.get_flux_total(band, True)

		return self.catalog_combined

		def get_flux_total(self, band, rest_frame=False, keys=["disk", "bulge", "point"]):

		# NOTE: override rest-frame behavior as EGG does not have absmags for disk/bulge separately
		if rest_frame:
		keys = ["", "point"]

		ret = np.zeros(self.catalog_combined.size)

		for p in keys:
		key = "magabs_" * rest_frame + f"{band}" + f"_{p}" * (p != "")
		flux = self.catalog_combined[key]
		if rest_frame:
		flux = util.mag_to_flux(flux)
		ret += np.where(np.isfinite(flux), flux, 0.0)

		if rest_frame:
		ret = util.flux_to_mag(ret)

		return ret

		def write(self):
		fitsio.write(self.filename, self.catalog_combined, clobber=True)
		return self.catalog_combined

		def get_is_galaxy(self):
		return np.isfinite(self.catalog_combined["Z"])

		def get_is_star(self):
		return ~self.get_is_galaxy()

		def get_index_star(self, i):

		"""
		Return the index in the stellar catalog corresponding to 'i' in the
		combined catalog
		"""

		n_galaxy, n_star, n_binary = self.get_number_galaxy_star_binary()

		# Handles single star
		if n_galaxy <= i < n_galaxy + n_star:
		return i - n_galaxy

		# Handles binary star
		if n_galaxy + n_star <= i:
		return i - n_galaxy - n_star

		# Handles something else
		return None


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

src/python/catalog_galaxy.py

0 → 100644

+150 −0

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-02-13 16:20:51

		"""
		Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution
		"""

		import logging
		import os

		import fitsio
		import numpy as np


		logger = logging.getLogger(__name__)

		class CatalogGalaxy:

		def __init__(self, dirname, egg):

		self.dirname = dirname
		self.egg = egg
		self.bands = [s.strip() for s in self.egg["BANDS"][0]]
		self.rfbands = [s.strip() for s in self.egg["RFBANDS"][0]]

		self.catalog = self.get_catalog()

		def get_catalog(self):

		# Create the Galaxy catalog
		self.catalog = np.empty_like(self.egg["RA"][0], dtype=self.get_dtype())
		filename = f"{self.dirname}/galaxy.fits"

		if not os.path.exists(filename):

		# column, type, description, unit
		for c, t, d, u in self.get_columns(self.bands, self.rfbands):
		self.catalog[c] = self.get_galaxy(c).astype(t)

		if not os.path.exists(os.path.dirname(filename)):
		os.makedirs(os.path.dirname(filename))

		fitsio.write(filename, self.catalog, clobber=True)

		self.catalog = fitsio.read(filename)
		return self.catalog

		@staticmethod
		def get_columns(bands, rfbands, is_public=False):

		ret = [

		# Columns from EGG
		("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"),
		("M", np.float64, "log10 of stellar mass", "Msun"),
		("SFR", np.float64, "Star-formation rate", "Msun/yr"),
		("PASSIVE", bool, "Is passive (non-star-forming)", ""),
		("CLUSTERED", bool, "Is clustered", ""),

		] + [

		# Galaxy morphological parameters
		("DISK_ANGLE", np.float64, "Galaxy disk rotation angle", "deg"),
		("DISK_RADIUS", np.float64, "Galaxy disk half-light radius", "arcsec"),
		("DISK_RATIO", np.float64, "Galaxy disk ratio of minor to major axis", ""),
		("BULGE_ANGLE", np.float64, "Galaxy bulge rotation angle", "deg"),
		("BULGE_RADIUS", np.float64, "Galaxy bulge half-light radius", "arcsec"),
		("BULGE_RATIO", np.float64, "Galaxy bulge ratio of minor to major axis", ""),

		] + [

		# Emission line extinction
		("AVLINES_BULGE", np.float64, "Emission line extinction in the bulge", "mag"),
		("AVLINES_DISK", np.float64, "Emission line extinction in the disk", "mag"),

		] + [

		# Disk fluxes
		(b.strip() + "_disk", np.float64, f"Galaxy disk {b} flux", "uJy") for b in bands

		] + [

		# Bulge fluxes
		(b.strip() + "_bulge", np.float64, f"Galaxy bulge {b} flux", "uJy") for b in bands

		] + [

		# EGG absolute magnitudes
		("magabs_" + b.strip(), np.float64, f"Galaxy {b} absolute magnitude", "ABmag") for b in rfbands

		]
		return ret

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

		def get_galaxy(self, key):

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

		_bands = [b.strip() for b in self.bands]
		key_without_suffix = key.replace("_disk", "").replace("_bulge", "")

		try:
		if key_without_suffix in _bands:
		# 20250514
		#idx = _bands.index(key_without_suffix)
		idx = self.bands.index(key_without_suffix)
		key_egg = "FLUX"
		if "bulge" in key:
		key_egg += "_BULGE"
		if "disk" in key:
		key_egg += "_DISK"
		return self.egg[key_egg][0, :, idx]

		elif "magabs_" in key:
		idx = self.rfbands.index(key.replace("magabs_", ""))
		return self.egg["RFMAG"][0, :, idx]
		except IndexError:
		# 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]
		if key == "M":
		# NOTE: Salpeter to Chabrier IMF conversion see e.g. Grylls+20
		ret -= 0.24

		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
		# of 0.63 (Chabrier) or 0.67 (Kroupa).

		# 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

		def __getitem__(self, key):
		return self.get_galaxy(key)

src/python/catalog_star.py

0 → 100644

+225 −0

File added.

Preview size limit exceeded, changes collapsed.