romanisim.models (#328)

* add romanisim.models

* in romanisim.models.wcs: change wcs.numerical_inverse() to wcs.invert()

* update import models.wcs, models.parameters

* modify romanisim.models APIs, now they also takes image_mod, reffiles as valid inputs

* modify texts in example_APIs notebook

* fix a wcs call error in image.py

* mv bandpass name conversions into models.bandpass. Fix a duplicated file read in models.bandpass

* fix a coner case in models.nonlinearity

* update models.bandpass to stick to SOC filter naming convention. Pass test_image, test_l1 unit tests

* update example API notebook

* add galsim_filter_name boolean option to read_gsfc_effarea in models.bandpass to make it also compatible with SOC naming convention.

* remove all galsim.roman dependence, replace the calls of romanisim.(parameters, wcs, psf, bandpass) with romanisim.models.(parameters, wcs, psf, bandpass)

* fix a  bug in models.wcs

* update integral nonlinearity correction to the current implementation in romanisim

* for reading files in the roman-tech-info, change:  RuntimeError to Exception. Add this try except block to models.bandpass

* update example notebook

* fix a bug in read_gsfc_effarea, remove unused import of asdf astropy.unit etc.

* Squashed 'romanisim/data/roman-technical-information/' content from commit 229bd3c

git-subtree-dir: romanisim/data/roman-technical-information
git-subtree-split: 229bd3ccc23e1a37fba76b62ed537a24d31f2e13

* add roman-technical-information repo as a subtree in romanisim.data, fix galsim.roman calls in docstring examples

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* ignore pre-commit check for roman-technical-information subtree

* fixed a typo in read_noise.py

* removed redundant fermi parameters

* default values for gain, dark_current, and read_noise are now retrieved from models.parameters, and can be updated via passing a metadata dict to the constructor

* moved romanisim.models.psf to romanisim.psf, modified the corresponding imports

* add __all__ to romanisim.psf

* import romanisim in bandpass docstring

* add a helper function to get the reference files (from local or crds)

* deprecated the contents in romanisim.wcs and romanisim.bandpass; corrct for an typo in the docstring in romanisim.models.bandpass

* update models.nonlinearity so it takes a integralnonlinearity datamodel as an input. Correspondingly updated the implementations in tests (test_l1 and test_linear).

* Added deprecation warnings in romanisim.nonlinearity. Also introduced backward compatibility in romanisim.models.nonlinearity for the functions repair_coefficients and evaluate_nl_polynomial. Note that romanisim.nonlinearity.NL is deprecated and is replaced by romanisim.models.nonlinearity.Nonlinearity.

* added a README.md for romanisim.data; changed the logic of accessing the default per SCA throughputs (now the defaults point to roman-tech-info local copy)

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* removed the unused throughput files in romanisim.data

* expose psf_utils and wcs_utils in romanisim.models.__init__

* fix a typo in nonlinearity, add reciprocity failure method in nonlinearity

* update bandpass to allow for getting a subset of Roman bandpasses

* add a getBandpass function (for getting a single band) to models.bandpass

* add relevant galsim.roman unittests

* fixed bugs in image.py and docstring

* fix style failures

* update docstrings for model classes

* remove romanisim/models/example_APIs.ipynb

---------

Co-authored-by: Eddie Schlafly <eschlafly@stsci.edu>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: 3 people

.gitignore

@@ -149,3 +149,4 @@ docs/source/api
 
 # misc
 uv.lock
+romanisim/models/example_APIs.ipynb

.pre-commit-config.yaml

@@ -1,7 +1,12 @@
 ci:
   autoupdate_schedule: monthly
 
-exclude: ".*\\.asdf$"
+# exclude: ".*\\.asdf$"
+exclude: |
+  (?x)^(
+      .*\\.asdf$|
+      romanisim/data/roman-technical-information/
+  )
 
 repos:
   # basic checks

romanisim/bandpass.py

@@ -5,211 +5,26 @@
 the nominal flat AB spectrum of 3631 Jy.  The ultimate source of this
 information is https://roman.gsfc.nasa.gov/science/WFI_technical.html .
 """
-import numpy as np
-from importlib import resources
-from scipy import integrate
-from astropy import constants
-from astropy.io import ascii
-from astropy import units as u
-from romanisim import parameters
 
-# to go from calibrated fluxes in maggies to counts in the Roman bands
-# we need to multiply by a constant determined by the AB magnitude
-# system and the shape of the Roman bands.
-# The constant is \int (3631 Jy) (1/hv) T(v) dv
-# T(v) should be the effective area at each wavelength, I guess
-# divided by some nominal overall effective area.
+import warnings
 
+warnings.warn(
+    "romanisim.bandpass is deprecated and will be removed in a future version. "
+    "Please use 'romanisim.models.bandpass' instead.",
+    DeprecationWarning,
+    stacklevel=2,
+)
 
-# provide some translation dictionaries for the mapping from
-# the galsim bandpass names to the Roman bandpass names and vice versa.
-# it would be nice to be agnostic about which one we use.
-galsim_bandpasses = [
-    'R062', 'Z087', 'Y106', 'J129', 'H158', 'F184', 'K213', 'W146']
-galsim2roman_bandpass = {x: 'F' + x[1:] for x in galsim_bandpasses}
-roman2galsim_bandpass = {v: k for k, v in galsim2roman_bandpass.items()}
+def __getattr__(name):
+    if name.startswith("__"):
+        raise AttributeError(name)
 
-# provide some no-ops if we are given a key in the right bandpass
-galsim2roman_bandpass.update(**{k: k for k in roman2galsim_bandpass})
-roman2galsim_bandpass.update(**{k: k for k in galsim_bandpasses})
+    from romanisim.models import bandpass
 
-# AB Zero Spectral Flux Density
-ABZeroSpFluxDens = 3631e-23 * u.erg / (u.s * u.cm**2 * u.hertz)
-
-
-def read_gsfc_effarea(sca, filename=None):
-    """Read an effective area file from Roman.
-
-    This just puts together the right invocation to get an Excel-converted
-    ECSV file into memory.
-
-    Parameters
-    ----------
-    sca : int
-        the name of the detector. A number between 1-18. 
-    filename : str
-        filename to read in
-
-    Returns
-    -------
-    astropy.table.Table
-        table with effective areas for different Roman bandpasses.
-    """
-
-    # The throughput files come as ECSV from March 2024 so the invocation changed to properly read the ECSVs. 
-    if filename is None:
-        filename = str(resources.files('romanisim') / 'data' / 'Roman_effarea_tables_20240327' 
-               / f'Roman_effarea_v8_SCA{sca:02}_20240301.ecsv')
-        out = ascii.read(filename)
-    else:
-        out = ascii.read(filename)
-    return out
-
-
-def compute_abflux(sca, effarea=None):
-    """Compute the AB zero point fluxes for each filter.
-
-    How many electrons would a zeroth magnitude AB star deposit in
-    Roman's detectors in a second?
-
-    Parameters
-    ----------
-    sca : int
-        the name of the detector. A number between 1-18. 
-    effarea : astropy.Table.table
-        Table from GSFC with effective areas for each filter.
-
-    Returns
-    -------
-    dict[str] : float
-        lookup table of zero point fluxes for each filter (electrons / s)
-    """
-
-    if effarea is None:
-        effarea = read_gsfc_effarea(sca)
-
-    # get the filter names.  This is a bit ugly since there's also
-    # a wavelength column 'Wave', and Excel appends a column to each line
-    # which astropy then gives a default name col12 to.
-    filter_names = [x for x in effarea.dtype.names
-                    if x != 'Wave' and 'col' not in x]
-    abfv = ABZeroSpFluxDens
-    out = dict()
-    for bandpass in filter_names:
-        out[bandpass] = compute_count_rate(flux=abfv, bandpass=bandpass, sca=sca, effarea=effarea)
-
-    # Saving the SCA information to use the correct throughput curves for each detector. 
-    out = {f'SCA{sca:02}':out}  
-    return out
-
-
-def get_abflux(bandpass, sca):
-    """Get the zero point flux for a particular bandpass.
-
-    This is a simple wrapper for compute_abflux, caching the results.
-
-    Parameters
-    ----------
-    bandpass : str
-        the name of the bandpass
-    sca : int
-        the name of the detector. A number between 1-18. 
-    Returns
-    -------
-    float
-        the zero point flux (electrons / s)
-    """
-    bandpass = galsim2roman_bandpass.get(bandpass, bandpass)
-
-    # If abflux for this bandpass for the specified SCA has been calculated, return the calculated
-    # value instead of rerunning an expensive calculation
-    abflux = getattr(get_abflux, 'abflux', None)
-    if (abflux is None) or (f'SCA{sca:02}' not in abflux):
-        abflux = compute_abflux(sca)
-        get_abflux.abflux = abflux
-    return abflux[f'SCA{sca:02}'][bandpass]
-
-
-def compute_count_rate(flux, bandpass, sca, filename=None, effarea=None, wavedist=None):
-    """Compute the count rate in a given filter, for a specified SED.
-
-    How many electrons would an object with SED given by
-    flux deposit in Roman's detectors in a second?
-
-    Parameters
-    ----------
-    flux : float or np.ndarray with shape matching wavedist.
-        Spectral flux density in units of ergs per second * hertz * cm^2
-    bandpass : str
-        the name of the bandpass
-    sca : int
-        the name of the detector. A number between 1-18. 
-    filename : str
-        filename to read in
-    effarea : astropy.Table.table
-        Table from GSFC with effective areas for each filter.
-    wavedist : numpy.ndarray
-        Array of wavelengths along which spectral flux densities are defined in microns
-
-    Returns
-    -------
-    float
-        the total bandpass flux (electrons / s)
-    """
-    # Read in default Roman effective areas from Goddard, if areas not supplied
-    if effarea is None:
-        effarea = read_gsfc_effarea(sca, filename)
-
-    # If wavelength distribution is supplied, interpolate flux and area
-    # over it and the effective area table layout
-    if wavedist is not None:
-        # Ensure that wavedist and flux have the same shape
-        if wavedist.shape != flux.shape:
-            raise ValueError('wavedist and flux must have identical shapes!')
-
-        all_wavel = np.unique(np.concatenate((effarea['Wave'], wavedist)))
-        all_flux = np.interp(all_wavel, wavedist, flux)
-        all_effarea = np.interp(all_wavel, effarea['Wave'], effarea[bandpass])
-    else:
-        all_wavel = effarea['Wave']
-        all_flux = flux
-        all_effarea = effarea[bandpass]
-
-    integrand = all_flux * constants.c / (
-        all_wavel * u.micron)**2  # f_lambda
-    integrand /= constants.h * constants.c / (
-        all_wavel * u.micron)  # hc/lambda
-    integrand *= all_effarea * u.m**2  # effective area in filter
-    # integrate.simpson looks like it loses units.  So convert to something
-    # we know about.
-    integrand = integrand.to(1 / (u.s * u.micron)).value
-
-    zpflux = integrate.simpson(integrand, x=all_wavel)
-    # effarea['Wave'] is in microns, so we're left with a number of counts
-    # per second
-
-    return zpflux
-
-
-def etomjysr(bandpass, sca):
-    """Compute factor converting e/s/pix to MJy/sr.
-
-    Assumes a pixel scale of 0.11" (romanisim.parameters.pixel_scale)
-
-    Parameters
-    ----------
-    bandpass : str
-        the name of the bandpass
-    sca : int
-        the name of the detector. A number between 1-18. 
-    Returns
-    -------
-    float
-        the factor F such that MJy / sr = F * DN/s
-    """
-
-    abflux = get_abflux(bandpass, sca)  # e/s corresponding to 3631 Jy
-    srpix = ((parameters.pixel_scale * u.arcsec) ** 2).to(u.sr).value
-    mjysr = 1 / abflux * 3631 / 10 ** 6 / srpix  # should be ~0.3
-    return mjysr
-    
+    warnings.warn(
+        f"'romanisim.bandpass.{name}' is deprecated and will be removed in a future version."
+        f"Please use 'romanisim.models.bandpass.{name}' instead.",
+        DeprecationWarning,
+        stacklevel=2,
+    )
+    return getattr(bandpass, name)

romanisim/catalog.py

@@ -7,16 +7,16 @@
 import dataclasses
 import numpy as np
 import galsim
-from galsim import roman
 from astropy import coordinates, table
 from astropy import units as u
 from astropy.io import fits
 import astropy_healpix
 import astropy.time
 from romanisim import gaia as rsim_gaia
-from . import util, log, parameters
-import romanisim.bandpass
+from . import util, log
+import romanisim.models.bandpass
 import yaml
+from romanisim.models import parameters
 
 # COSMOS constants taken from the COSMOS2020 paper:
 # https://arxiv.org/pdf/2110.13923
@@ -32,7 +32,7 @@
 F184_KS_COEFF = 0.3838145747397368
 
 # Bandpass filters
-BANDPASSES = list(romanisim.bandpass.galsim2roman_bandpass.values())
+BANDPASSES = set(romanisim.models.bandpass.galsim2roman_bandpass.values())
 
 
 @dataclasses.dataclass
@@ -83,20 +83,20 @@ def make_dummy_catalog(coord,
         rng = galsim.UniformDeviate(seed)
 
     if galaxy_sample_file_name is None:
-        cat1 = galsim.COSMOSCatalog(sample='25.2', area=roman.collecting_area,
+        cat1 = galsim.COSMOSCatalog(sample='25.2', area=parameters.collecting_area,
                                     exptime=1)
-        cat2 = galsim.COSMOSCatalog(sample='23.5', area=roman.collecting_area,
+        cat2 = galsim.COSMOSCatalog(sample='23.5', area=parameters.collecting_area,
                                     exptime=1)
     else:
         cat1 = galsim.COSMOSCatalog(galaxy_sample_file_name,
-                                    area=roman.collecting_area, exptime=1)
+                                    area=parameters.collecting_area, exptime=1)
         cat2 = cat1
 
     if chromatic:
         # following Roman demo13, all stars currently have the SED of Vega.
         # fluxes are set to have a specific value in the y bandpass.
         vega_sed = galsim.SED('vega.txt', 'nm', 'flambda')
-        y_bandpass = roman.getBandpasses(AB_zeropoint=True)['Y106']
+        y_bandpass = romanisim.models.bandpass.getBandpasses(AB_zeropoint=True)['Y106']
 
     objlist = []
     locs = util.random_points_in_cap(coord, radius, nobj, rng=rng)
@@ -115,7 +115,7 @@ def make_dummy_catalog(coord,
             mu = np.log(mu_x**2 / (mu_x**2 + sigma_x**2)**0.5)
             sigma = (np.log(1 + sigma_x**2 / mu_x**2))**0.5
             gd = galsim.GaussianDeviate(rng, mean=mu, sigma=sigma)
-            flux = np.exp(gd()) / roman.exptime
+            flux = np.exp(gd()) / parameters.exptime
             if chromatic:
                 sed = vega_sed.withFlux(flux, y_bandpass)
                 obj = galsim.DeltaFunction() * sed
@@ -411,8 +411,8 @@ def make_galaxies(coord,
         Table for use with table_to_catalog to generate catalog for simulation.
     """
     if bandpasses is None:
-        bandpasses = roman.getBandpasses().keys()
-        bandpasses = [romanisim.bandpass.galsim2roman_bandpass[b]
+        bandpasses = romanisim.models.bandpass.getBandpasses().keys()
+        bandpasses = [romanisim.models.bandpass.galsim2roman_bandpass[b]
                       for b in bandpasses]
     if rng is None:
         rng = galsim.UniformDeviate(seed)
@@ -592,8 +592,8 @@ def make_stars(coord,
         Table for use with table_to_catalog to generate catalog for simulation.
     """
     if bandpasses is None:
-        bandpasses = roman.getBandpasses().keys()
-        bandpasses = [romanisim.bandpass.galsim2roman_bandpass[b]
+        bandpasses = romanisim.models.bandpass.getBandpasses().keys()
+        bandpasses = [romanisim.models.bandpass.galsim2roman_bandpass[b]
                       for b in bandpasses]
     if rng is None:
         rng = galsim.UniformDeviate(seed)

romanisim/cr.py

@@ -1,7 +1,7 @@
 import numpy as np
 import scipy.interpolate as interpolate
 
-from . import parameters
+from romanisim.models import parameters
 
 
 def create_sampler(pdf, x):