Add utility and demonstration scripts

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2025-07-03 19:33:15

"""Test the combined catalog."""

import logging

import os

import numpy as np

from lsst_inaf_agile.catalog_agn import CatalogAGN

from lsst_inaf_agile.catalog_combined import CatalogCombined

from lsst_inaf_agile.catalog_galaxy import CatalogGalaxy

from lsst_inaf_agile.catalog_star import CatalogStar

from lsst_inaf_agile.egg import Egg

from lsst_inaf_agile.image_simulator import ImageSimulator

from lsst_inaf_agile.merloni2014 import Merloni2014

logger = logging.getLogger(__name__)

def get_catalog_combined(dirname):

    """Return a test CatalogCombined object."""

    os.makedirs(dirname, exist_ok=True)

    # Initialize EGG

    egg = Egg.run_config("etc/config_egg.ini")

    # Create the individual catalogs

    catalog_galaxy = CatalogGalaxy(dirname, egg)

    catalog_agn = CatalogAGN(

        dirname=dirname,

        catalog_galaxy=catalog_galaxy,

        type_plambda="zou+2024",

        save_sed=1,

        seed=20250621,

        merloni2014=Merloni2014(1, 0, 0.05, 0.95),

        filter_db="data/egg/share/filter-db/db.dat",

    )

    catalog_star = CatalogStar(dirname, catalog_galaxy, is_binary=False)

    catalog_binary = CatalogStar(dirname, catalog_galaxy, is_binary=True)

    # Create the combined catalog

    catalog_combined = CatalogCombined(dirname, catalog_galaxy, catalog_agn, catalog_star, catalog_binary)

    assert np.isclose(catalog_combined.get_area(), 0.01)

    return catalog_combined

def main():

    # Setup the dirname

    dirname = "data/catalog/test"

    # Get the combined catalog

    catalog_combined = get_catalog_combined(dirname)

    # Create the reference catalog

    maglim = 24

    selection_band = "lsst-r"

    catalog_combined.write_reference_catalog(f"{dirname}/reference_catalog.csv", maglim, selection_band)

    # Simulate images

    image_simulator = ImageSimulator(

        f"{dirname}/imsim", catalog_combined, "data/baseline/baseline_v4.0_10yrs.db"

    )

    # NOTE: select the first visit from COSMOS DDF

    visits = image_simulator.get_visit()

    select = visits["target_name"] == "DD:COSMOS"

    observation_id = visits["observationId"][select].iloc[0]

    ra_dec = +150.11916667, +2.20583333

    image_simulator.write_instance_catalog(observation_id, ra_dec=ra_dec)

    image_simulator.simulate_image(observation_id, detector=94)

if __name__ == "__main__":

    main()

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2025-10-24 20:30:23

"""

Create a catalog of lightcurve outliers for interactive analysis

"""

import os

import sqlite3

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd

from astropy.table import Table

if not os.path.exists("catalog.fits"):

    with sqlite3.connect("data/catalog/dr1_new_new/db/20250725/master.db") as con:

        df = pd.read_sql_query(

            """

            SELECT

                f.objectId,

                f.ccdVisitId,

                o.r_ra,

                o.r_dec,

                o.r_psfFlux,

                psfFlux,

                psfFluxErr,

                zenithDistance,

                nPsfStar

            FROM ForcedSource AS f

            JOIN Object AS o USING (objectId)

            JOIN CcdVisit USING (ccdVisitId)

            WHERE

                o.detect_isPrimary = 1 AND

                o.r_psfFlux > 912 AND

                o.r_psfFlux_flag = 0 AND

                o.r_extendedness_flag = 0 AND

                o.r_extendedness = 0 AND

                f.psfFlux_flag = 0 AND

                f.band = 'r'

            """,

            con,

        )

    from astropy.table import Table

    table = Table.from_pandas(df)

    object_ids = np.unique(table["objectId"])

    is_outlier = np.empty_like(table["objectId"], dtype=bool)

    for i, object_id in enumerate(object_ids):

        print(i, object_ids.size)

        select = table["objectId"] == object_id

        rows = table[select]

        fluxes = table["psfFlux"][select]

        flux = table["r_psfFlux"][select]

        median = np.median(fluxes)

        std = np.median(np.abs(fluxes - median))

        is_outlier[select] = (select.sum() > 30) * np.abs(fluxes - median) > 5 * std

    table["is_outlier"] = is_outlier

    table.write("catalog.fits", overwrite=True)

from astropy.table import Table

table = Table.read("catalog.fits")

bins = np.arange(0, 90, 10)

cens = (bins[:-1] + bins[1:]) * 0.5

fig, axes = plt.subplots(2, 1)

hist1 = np.histogram(table["zenithDistance"], bins=bins)[0]

hist2 = np.histogram(table["zenithDistance"][table["is_outlier"]], bins=bins)[0]

ratio = np.ma.true_divide(hist2, hist1)

axes[0].plot(cens, ratio)

bins = np.arange(0, 1000, 10)

cens = (bins[:-1] + bins[1:]) * 0.5

hist1 = np.histogram(table["nPsfStar"], bins=bins)[0]

hist2 = np.histogram(table["nPsfStar"][table["is_outlier"]], bins=bins)[0]

ratio = np.ma.true_divide(hist2, hist1)

axes[1].plot(cens, ratio)

plt.show()

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2025-12-02 14:59:18

"""

Create IY AGILE dataset

"""

# Read in IDs

import os

import sqlite3

import numpy as np

import pandas as pd

from astropy.table import Table

ids: list[int] = []

nrandom = 100

np.random.seed(20251202)

for _, where in [

    ("AGN1", "Z > 0 AND is_agn = 1 AND is_optical_type2 = 0"),

    ("AGN2", "Z > 0 AND is_agn = 1 AND is_optical_type2 = 1"),

    ("Galaxy", "Z > 0 AND is_agn = 0"),

]:

    with sqlite3.connect("/media/aetviita/Uranus/data/agile/dr1/db/20250725/master.db") as con:

        df = pd.read_sql_query(

            f"""

            SELECT ID

            FROM Truth

            WHERE {where}

            LIMIT 100

            """,

            con,

        )

        t = Table.from_pandas(df)

    print("\n".join(str(i) for i in t["ID"]))

    continue

    for ID in t["ID"]:

        dirname = "data/example_seds"

        if not os.path.exists(dirname):

            os.makedirs(dirname)

        for component in "bulge", "disk":

            cmd = (

                f"egg-getsed component={component}"

                " seds=/media/aetviita/Uranus/data/agile/dr1/egg-seds.dat id={ID}"

                " out={dirname}/egg-seds-{component}-{ID}.fits"

            )

            os.system(cmd)

        filename = f"/media/aetviita/Uranus/data/agile/dr1/seds/agn-seds-{ID}.fits"

        if os.path.exists(filename):

            os.system(f"cp -v {filename} data/example_seds/")

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2025-11-28 15:19:18

"""

Draw a lightcurve animation

"""

import fitsio

import matplotlib.pyplot as plt

import numpy as np

from matplotlib.animation import FuncAnimation

fig, axes = plt.subplots(1, 2)

xdata, ydata = [], []

images = [fitsio.read("notebooks/cutout_calexp.fits") for i in range(10)]

img = axes[0].imshow(images[0])

(ln,) = axes[1].plot([], [], "ro")

def init():

    axes[1].set_xlim(0, 2 * np.pi)

    axes[1].set_ylim(-1, 1)

    return (ln,)

def update(frame):

    img.set_data(images[frame])

    xdata.append(frame / 100 * 2 * np.pi)

    ydata.append(np.sin(frame / 100 * 2 * np.pi))

    ln.set_data(xdata, ydata)

    return ln, img

ani = FuncAnimation(fig, update, frames=np.arange(100), init_func=init, blit=True)

plt.show()