Minor updates

[common]

# The simulated survey area in square degrees. Note that the area is used

# to set the directory name for many of the outputs. Look for ${area} entries

#save_sed: 0

save_sed: 1

area: 0.01

#area: 0.04

#area: 0.10

#area: 0.50

#area: 1.00

#area: 2.00

#area: 12.00

#area: 15.00

#area: 18.00

#area: 24.00

#area: 100.00

zmin: 0.21

zmax: 5.49

dirname_prefix: data/catalog/test.frozen

log: ${common:dirname_prefix}/run.log

seds_file: ${common:dirname_prefix}/seds

bands:   [lsst-u,johnson-B,lsst-g,lsst-r,lsst-i,lsst-z,lsst-y,euclid-VIS,euclid-nisp-Y,euclid-nisp-J,euclid-nisp-H,spitzer-irac1,spitzer-irac2,spitzer-irac3,spitzer-irac4,wise-w1,wise-w2,wise-w3,wise-w4,mock-1000-4000,mock-1450,mock-4400,mock-5007,mock-15um,vista-Y,vista-J,vista-H,vista-Ks]

rfbands: ${common:bands}

[egg]

# EGG filename

out: ${common:dirname_prefix}/egg.fits

# NOTE: these arguments are fed verbatim to the egg. See egg-gencat help for

# more information

area: ${common:area}

verbose: 1

#mass_func: opt/egg/share/mass_func_cosmos2020_agn_zou24_loglambda32.fits

#mass_func: opt/egg/share/mass_func_candels.fits

mass_func: data/egg/share/mass_func_cosmos2020_agn_zou24_loglambda32.fits

## NOTE: COSMOS r-band 10yr limiting magnitude is 28.55 (baseline v3.4, http://astro-lsst-01.astro.washington.edu:8080/summaryStats?runId=4)

## Simulating 1.5mag deeper should be sufficient

#maglim: 30

#selection_band: lsst-r

zmin: ${common:zmin}

zmax: ${common:zmax}

# NOTE: The stellar mass input argument is in Salpter IMF

#mmin: 8.00

#mmin: 8.50

mmin: 8.74

#mmin: 9.00

#mmin: 9.24

#mmin: 9.50

#mmin: 10.50

dz: 0.05

# NOTE: 0.30 deg should cover most cases i.e. r~0-200 Mpc/h at z~0.20-5.50

clust_r0: 0.30

clust_fclust_him: 0.40

clust_fclust_lom: 0.20

bands: ${common:bands}

rfbands: ${common:bands}

ra0: +150.11916667

dec0: +2.20583333

save_sed: ${common:save_sed}

seed: 20230412

no_flux: 0

filter_db: data/egg/share/filter-db/db.dat

[combined]

# Area simulated by EGG

area: ${common:area}

# Band simulated by egg

bands: ${common:bands}

rfbands: ${common:bands}

# Filename of the combined output catalog

filename_catalog: ${common:dirname_prefix}/catalog.fits

# Filename of the AGN catalog

filename_agn: ${common:dirname_prefix}/agn.fits

# Filename of the EGG catalog

filename_egg: ${common:dirname_prefix}/egg.fits

# Filename of the stars catalog

filename_stars: ${common:dirname_prefix}/stars.fits

# Filename of the binary star catalog

filename_binaries: ${common:dirname_prefix}/binaries.fits

# Directory to save the SEDs in

seds_file: ${common:seds_file}

# Plambda to use

#type_plambda: bongiorno+2016

type_plambda: zou+2024

# Whether to use classification or no

no_classification: 1

# delta mass for grid in stellar mass

dm: 0.10

# delta redshift for grid in redshift

dz: 0.10

# delta LX for grid in luminosity

dl: 0.10

# Save AGN SEDs or no

save_sed: ${common:save_sed}

## Random seed

#seed: 20240326

# Fetch stars or no

fetch_stars: 1

#fetch_stars: 0

# Parameters for Merloni+2014 optical type assignment

merloni2014_interpolate: 0

merloni2014_extrapolate: 0

merloni2014_f_obs_minimum: 0.00

merloni2014_f_obs_maximum: 1.00

[reference_catalog]

# Filename of the reference catalog

filename: ${common:dirname_prefix}/reference_catalog.csv

# Magnitude limit for the reference catalog

#maglim: 20

maglim: 24

# Selection band for the magnitude limit

selection_band: lsst-r

[imsim]

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

# SQL query for the cadence database to define observationIds for running

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

# Updated 20240716: query can be a filename with one observationId per row

#query: data/observation_ids_cosmos_3yr.txt

#query: data/observation_ids_cosmos_3yr_gri.txt

#query: data/observation_ids_cosmos_3yr_last.txt

#query: data/observation_ids_cosmos_test.txt

#query: data/observation_ids_cosmos_3yr_gri_head10.txt

#query: data/observation_ids_cosmos_3yr_missing.txt

query: data/baseline_v4.0_observation_ids_cosmos_3yr.txt

# NOTE: all COSMOS

#query: SELECT * FROM observations WHERE note="DD:COSMOS"

## NOTE: r-band observations

#query: SELECT * FROM observations WHERE (filter='r' AND note="DD:COSMOS" AND observationId=29551)

##query: SELECT * FROM observations WHERE (filter='r' AND note="DD:COSMOS")

## NOTE: first 20 COSMOS r-band observations

#query: SELECT * FROM observations WHERE ((observationId>=29551 AND observationID<29554) OR (observationId>=30175 AND observationId<30192))

## NOTE: 20point r-band cadence

#query: SELECT * FROM observations WHERE (observationId=29551 OR observationId=30175 OR observationId=44921 OR observationId=45407 OR observationId=59084 OR observationId=72699 OR observationId=88519 OR observationId=103589 OR observationId=118569 OR observationId=145140 OR observationId=157286 OR observationId=234484 OR observationId=235890 OR observationId=236499 OR observationId=241887 OR observationId=242456 OR observationId=243196 OR observationId=243272 OR observationId=243862 OR observationId=243950)

## NOTE: first observation in each band. The order of the first sequence is giruyz

#query: SELECT * FROM observations WHERE (observationId=29521 OR observationId=29531 OR observationId=29551 OR observationId=30192 OR observationId=30200 OR observationId=59122)

## NOTE: first observation of gri bands

#query: SELECT * FROM observations WHERE (observationId=29521 OR observationId=29531 OR observationId=29551)

## NOTE: first six of each sequence

#query: SELECT * FROM observations WHERE (observationId=30192 OR observationId=45409 OR observationId=88539 OR observationId=103609 OR observationId=157306 OR observationId=242461 OR observationId=29521 OR observationId=44891 OR observationId=59054 OR observationId=72669 OR observationId=88489 OR observationId=103559 OR observationId=29551 OR observationId=44921 OR observationId=59084 OR observationId=72699 OR observationId=88519 OR observationId=103589 OR observationId=29531 OR observationId=44901 OR observationId=59064 OR observationId=72679 OR observationId=88499 OR observationId=103569 OR observationId=59122 OR observationId=72737 OR observationId=118607 OR observationId=145178 OR observationId=235928 OR observationId=236537 OR observationId=30200 OR observationId=45417 OR observationId=59104 OR observationId=72719 OR observationId=88547 OR observationId=103617)

# Magnitude limit for the instance catalog

maglim: 30

#maglim: 24

# Selection band for the magnitude limit

selection_band: lsst-r

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

# List of detectors to simulate. Number ordering follows the scheme here:

# 	https://lsstdesc.org/imSim/lsst-camera.html#specifying-which-sensors-to-simulate

# 	http://lsstdesc.org/imSim/lsst-camera.html

# 1det center

#det_list: 94

# 9det center

#

#	096 097 098

#	093 094 095

#	090 091 092

#

#det_list: 90 91 92 93 94 95 96 97 98

# 21det center leaving out the corners:

#

# 	    135 136 137

# 	089 096 097 098 105

# 	086 093 094 095 102

# 	083 090 091 092 099

# 	    051 052 053

#

det_list: 51 52 53 83 86 89 90 91 92 93 94 95 96 97 98 99 102 105 135 136 137

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

# Based directory name for the imsim output products

dirname: ${common:dirname_prefix}/imsim

#dirname: ${common:dirname_prefix}/imsim_gri

#dirname: ${common:dirname_prefix}/imsim_test_reff

# Visits definition file

#filename_visits: /staff2/viitanen/agile/data/baseline_v3.0_10yrs.db

#filename_visits: /staff2/viitanen/agile/data/baseline_v3.0_10yrs.db

# White-space separated list of exptimes for ugrizy, in sequence. 30s is the standard single visit

# exposure time. The second option simulates the whole single visit DDF

# sequence (ugrizy): (8, 10, 20, 20, 24, 18) exposures.

exptime: 30 30 30 30 30 30

#exptime: 240 300 600 600 720 540

# Whether to estimate AGN lightcurves or not

do_lightcurve: 1

#do_lightcurve: 0

# Whether to run galsim

#run_galsim: 1

run_galsim: 0

# Turn on AGN host galaxies?

write_agn_host_galaxy: 1

#write_agn_host_galaxy: 0

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

# NOTE: the following will override default options from the simulation

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

# Override ra pointing

ra: 150.11916667

# Override dec pointing

dec: 2.20583333

[figure]

dirname: ${common:dirname_prefix}/figure

  - astropy

  - batman-package

  - fitsio

  - jupyter

  - matplotlib

  - numpy

  - scipy

from copy import deepcopy

from itertools import product

import glob

from multiprocessing import Pool

import multiprocessing

import os

import sys

import time

def _get_lightcurve_agn(args):

    i, filename, band, flux, z, mag_i, logMbh, type2 = args

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

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

    kwargs = {

        "mjd0": 0,

                for i in np.arange(self.catalog.size)

            ]

            import multiprocessing

            cpu_count = multiprocessing.cpu_count()

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

            with Pool(processes=cpu_count) as p:

            t0 = time.time()

            with multiprocessing.Pool(processes=cpu_count) as p:

                log_lambda_SAR = p.starmap(util.get_log_lambda_SAR, arguments)

            logger.info("... done")

            t1 = time.time()

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

        return log_lambda_SAR

                for z in np.arange(np.floor(self.catalog["Z"].min()), np.ceil(self.catalog["Z"].max()), self.dz):

                    for m in np.arange(np.floor(self.catalog["M"].min()), np.ceil(self.catalog["M"].max()), self.dm):

                        logger.info(f"{t} {z} {m}")

                        logger.debug(f"{t} {z} {m}")

                        select = self._get_select(z=z, m=m, t=t)

                        if select.sum() == 0:

                            continue

                                np.save(filename, func(10 ** (m + self.dm / 2), z + self.dz / 2, type=t))

                            ret += [np.load(filename)]

                        duty_cycle[select] = ret[1] / ret[0]

                        logger.info(f"Duty cycle t={t} z={z:6.2f} m={m:6.2f} N={select.sum()} f={ret[1]/ret[0]}")

                        logger.debug(f"Duty cycle t={t} z={z:6.2f} m={m:6.2f} N={select.sum()} f={ret[1]/ret[0]}")

        elif self.type_plambda == "zou+2024":

            duty_cycle = self.catalog["log_lambda_SAR"] >= 32.0

                    #   N_CTK = X * N_CTN / (1 - X)

                    select = self._get_select(z=z, l=l)

                    logger.info(f"{z} {l}")

                    logger.debug(f"{z} {l}")

                    if select.sum() == 0:

                        continue

                    frac = [ueda2014.get_f(l + self.dl / 2, z + self.dz / 2, nh) for nh in [20, 21, 22, 23, 24, 25]]

                    norm = np.sum(frac)

                    frac_agn_ctk = (frac[-2] + frac[-1]) / norm

                    logger.info(frac_agn_ctk)

                    logger.debug(frac_agn_ctk)

                    n_gal = select.sum()

            for i in range(self.catalog.size):

                logger.info(f"{i}  {self.catalog.size}  {100 * i / self.catalog.size}")

                logger.debug(f"{i}  {self.catalog.size}  {100 * i / self.catalog.size}")

                if self.catalog["log_lambda_SAR"][i] >= 32.0:

                    continue

    def _get_agn_sed(self, i, ratio_max=0.90):

        logger.info(f"Getting AGN SED {i} {self.catalog.size}")

        logger.debug(f"Getting AGN SED {i} {self.catalog.size}")

        # Get the wavlen in angstrom

        # NOTE: Add some interesting wavelengths for greater accuracy

            "seed": i + 1

        })

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

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

        import lightcurve

        lc = lightcurve.get_lightcurve_agn(*args, **kwargs)

            idxs = uid_agn

        idxs = np.atleast_1d(idxs)

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

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

        for ndone, i in enumerate(idxs):

            logger.info(f"{i} {ndone} {idxs.size} {100 * ndone / idxs.size:7.2}")

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

        filename = f"{self.dirname}/lightcurves/{i}/{band}.npy"

        if os.path.exists(filename):

            tt, lc = np.load(filename).T

            logger.info(f"Found binary star lightcurve {i} in {band}")

            logger.debug(f"Found binary star lightcurve {i} in {band}")

            return tt, lc

        # Try to create the lightcurve

                uid = c["ID"][i]

                idx = c["ID"][i] - (self.catalog.size - self.star_binary.size)

                s = self.star_binary[idx]

                logger.info(f"{uid} {idx} {i} {my_i} {select_binary.sum()} {100 * my_i / select_binary.sum()}")

                logger.debug(f"{uid} {idx} {i} {my_i} {select_binary.sum()} {100 * my_i / select_binary.sum()}")

                # NOTE: find out which one is the "primary" based on the radius

                i1, i2 = 1, 2

            f.write('\n')

        # Copy over the SED template

        os.system(f"cp -vn data/imsim/Constnu.Template.spec.gz {dirname_id}")

        os.system(f"cp -v --update=none data/imsim/Constnu.Template.spec.gz {dirname_id}")

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

        # Run galsim

Wrapper for calling egg

"""

import logging

logger = logging.getLogger()

import os

import subprocess

from util import ROOT

# See if EGG has been installed properly

        if not os.path.exists(dirname):

            os.makedirs(dirname)

        print("Running...", cmd)

        return os.system(cmd)

        logger.info(f"Running {cmd}")

        process = subprocess.Popen(

            cmd.split(),

            stdout=subprocess.PIPE,

            stderr=subprocess.PIPE,

            text=True,

            universal_newlines=True

        )

        for line in process.stdout:

            logger.info(f"{line.strip()}")

        for line in process.stderr:

            logger.info(f"{line.strip()}")

        return_code = process.wait()

        return return_code

    def get_sed(self, i):

        """Return a EGG SED"""

    # Simulate images

    # NOTE: requires imSim

    if False:

    if True:

        print("Simulating the images...")

        catalog.simulate_images(**config["imsim"])

    # Setup the logger

    logging.basicConfig(

        level=logging.INFO,

        format='%(asctime)s - %(name)s - %(levelname)s - %(message)s',

        #format='%(asctime)s - %(name)s - %(levelname)s - %(message)s',

        handlers=[

            logging.StreamHandler(),

            logging.FileHandler(config["common"]["log"])