Unverified Commit 2dd405fd authored by Akke Viitanen's avatar Akke Viitanen
Browse files

raise catalog_star coverage

parent 985e8ace

.pre-commit-config.yaml

+1 −0
Original line number Diff line number Diff line 
fail_fast: true



repos:

    # Compare the local template version to the latest remote template version

Makefile

+2 −4
Original line number Diff line number Diff line 
.PHONY: test

.PHONY: test tests



dr1:

	python3 python/main.py --config etc/config_dr1.ini
@@ -266,9 +266,7 @@ test: 
		python3 src/python/main.py --config etc/config_test.ini



tests:

	PYTHONPATH=opt:src/python:$(PYTHONPATH)

		python3 src/tests/test_truth_catalog.py



	python -m pytest --cov=./src --cov-report=html tests/



debug:

	PYTHONPATH=opt:src/python:$(PYTHONPATH) \

src/lsst_inaf_agile/catalog_star.py

+28 −10
Original line number Diff line number Diff line
@@ -12,6 +12,7 @@ import numpy as np 
import pyvo as vo

from astropy import constants

from astropy import units as u

from astropy.table import Table



from . import lightcurve, util
@@ -90,7 +91,7 @@ class CatalogStar: 
        """Return column corresponding to 'k'."""

        return self.catalog[k]



    def get_stars(self, selection="rmag", maglim=28, fbin=1.0):

    def get_stars(self, selection="rmag", maglim=28, debug=False):

        """Query NOIRlab to retrieve the stellar catalogs from LSST-SIM."""

        if os.path.exists(self.filename):

            return util.read_fits(self.filename)
@@ -114,17 +115,27 @@ class CatalogStar: 
                """

        logger.info(query)



        # Run the syncronous job

        tap_service = vo.dal.TAPService("https://datalab.noirlab.edu/tap")

        tap_results = tap_service.search(query)

        table = tap_results.to_table()



        import warnings



        from astropy.units import UnitsWarning



        with warnings.catch_warnings():

            warnings.simplefilter("ignore", UnitsWarning)



            if debug:

                # Return dummy output

                if self.is_binary:

                    table = Table.read("data/tests/test_catalog_star/bak/binaries.fits")

                else:

                    table = Table.read("data/tests/test_catalog_star/bak/stars.fits")

                table = table[table[selection] < maglim]



            else:

                # Run the syncronous job

                tap_service = vo.dal.TAPService("https://datalab.noirlab.edu/tap")

                tap_results = tap_service.search(query)

                table = tap_results.to_table()



            table.write(self.filename)



        return table
@@ -135,18 +146,24 @@ class CatalogStar: 
        pmode = self.stars["pmode"][i]

        return (label >= 4) & (label <= 6) & (pmode >= 0) & (pmode <= 1)



    def is_lpv(self, i, orich=True, crich=True, c_o=None):

    def is_c_rich(self, i=None):

        ret = self.stars["c_o"] > 1

        if i is None:

            return ret

        return ret[i]



    def is_lpv(self, i, orich=True, crich=True):

        """Return true if star 'i' is a long-period variable star."""

        label = self.stars["label"][i]

        pmode = self.stars["pmode"][i]

        select = (label >= 7) & (label <= 8) & (pmode >= 0) & (pmode <= 4)



        # Select crich/orich only?

        if crich ^ orich and c_o is not None:

        if crich ^ orich:

            if orich:

                select &= c_o <= 1

                select &= ~self.is_c_rich(i)

            if crich:

                select &= c_o > 1

                select &= self.is_c_rich(i)



        return select
@@ -201,6 +218,7 @@ class CatalogStar: 
        kwargs = {

            "mjd0": mjd0,

            "mjd": mjd,

            "band": band,

            "flux": self.catalog[f"{band}_point"][i],

            "pmode": pmode,

            "period": period,

src/lsst_inaf_agile/lightcurve.py

+21 −14
Original line number Diff line number Diff line
@@ -86,6 +86,9 @@ def get_lightcurve_agn( 


def get_lightcurve_cepheid(mjd0, mjd, flux, band, z, y, logte, logl, mass, pmode, period, seed=None):

    """Return a cepheid light-curve."""



    from CEP_lsst_quasar_project.get_lightcurve import get_lightcurve_normalized_cepheid



    source = {

        "z": z,

        "y": y,
@@ -97,16 +100,14 @@ def get_lightcurve_cepheid(mjd0, mjd, flux, band, z, y, logte, logl, mass, pmode 
        f"{band[-1]}mag": util.flux_to_mag(flux),

    }



    from CEP_lsst_quasar_project.get_lightcurve import get_lightcurve_normalized_cepheid



    time, flux = get_lightcurve_normalized_cepheid(source, band)



    # NOTE: add a random phase to the lightcurve

    np.random.seed(seed)

    time += np.random.uniform(0, period)

    time -= time.min()



    return np.interp(mjd - mjd0, time, flux, period=period)

    ret = np.interp((mjd - mjd0) % period, time, flux)

    return ret





def get_lightcurve_lpv(mjd0, mjd, flux, band, pmode, period, is_c_rich, seed=None):
@@ -142,9 +143,11 @@ def get_lightcurve_lpv(mjd0, mjd, flux, band, pmode, period, is_c_rich, seed=Non 
    tt = np.arange(0, period, 0.10) + np.random.uniform(0, period)

    A = get_amplitude_lpv()

    B = util.flux_to_mag(flux)

    yy = A * np.sin(2 * np.pi * tt / period) + B

    phi = (2 * np.pi * tt / period) % (2 * np.pi)

    yy = A * np.sin(phi) + B

    yy = util.mag_to_flux(yy)

    return np.interp(mjd - mjd0, tt, yy, period=period)



    return np.interp(mjd - mjd0, tt, yy, period=tt.max())





try:
@@ -156,6 +159,9 @@ except ImportError: 
    )





LCB = {}





def get_lightcurve_binary(

    mjd0,

    mjd,
@@ -214,17 +220,19 @@ def get_lightcurve_binary( 


    Returns

    -------

    lc1, lc2:

        lightcurves of the primary and the secondary. The total lightcurve of

        the system is then lc1 + lc2. In the same units of flux1 and flux2.

    lc: float

        Combined lightucrve of the primary and secondary in the same units as

        flux1 and flux2.



    """

    # NOTE: short-circuit on existing filename optimization

    if filename is not None and os.path.exists(filename) and not overwrite:

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

        if filename not in LCB:

            LCB[filename] = np.load(filename).T

        t, lc = LCB[filename]

        if return_full:

            return t, lc

        return np.interp(mjd, t, lc, period=True)

        return np.interp((mjd - mjd0) % p + t.min(), t, lc)



    # NOTE: there is some weird feature in batman with radius1 == radius2 which

    # is prevented by adding a very minor offset to the secondary radius
@@ -264,7 +272,7 @@ def get_lightcurve_binary( 
            return None



    # NOTE: debugging

    if False:

    if True:

        print("", radius1, flux1)

        print("", radius2, flux2)

        print("", param.t0, "time of inferior conjunction")
@@ -327,5 +335,4 @@ def get_lightcurve_binary( 


    # Interpolate the lightcurves

    # TODO: average over 30s?

    flux = np.interp(mjd, t, lc, period=True)

    return flux
 
    return np.interp((mjd - mjd0) % p + t.min(), t, lc)

src/lsst_inaf_agile/util.py

+8 −0
Original line number Diff line number Diff line
@@ -34,6 +34,14 @@ def mag_to_flux(mag): 
    return 3631 * 1e6 * 10 ** (mag / -2.5)





def mag_sum(mag):

    """

    Sum together magnitudes 'mag', return the combined magnitude.

    """

    flux = np.sum(mag_to_flux(mag))

    return flux_to_mag(flux)





def get_volume(zmin, zmax, area_deg2, H0=70.0, Om0=0.30, Tcmb0=2.73):

    """Return the comoving volume in Mpc for a redshift shell."""

    cosmo = FlatLambdaCDM(H0=H0, Om0=Om0, Tcmb0=Tcmb0)