Performance PR (#123)

* improved polynomial performance, fixing minor issues, standardized the fetching of meta data, and introducing numba

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

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

* added missing files, checked something for RAvila

* fixed doc string

---------

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

Author: Russell-Ryan

pyproject.toml

@@ -38,6 +38,7 @@ dependencies = [
     'h5py>=3.11',
     'drizzlepac>=3.9.1',
     'networkx>=3.3',
+    'numba>=0.62.0',
     'numpy>=2.0',
     'matplotlib>=3.9',
     'pandas>=2.2',

slitlessutils/core/modules/extract/multi/matrix.py

@@ -52,8 +52,12 @@ class Matrix:
     INT = np.uint64     # DO NOT CHANGE THIS
     FLOAT = np.float64  # DO NOT CHANGE THIS
 
-    def __init__(self, extorders, mskorders=None, invmethod='lsqr', path='tables',
-                 minunc=1e-10):
+    def __init__(self, extorders, **kwargs):
+        # extract the defaults
+        mskorders = kwargs.get('mskorders', None)
+        invmethod = kwargs.get('invmethod', 'lsqr')
+        path = kwargs.get('path', 'su_tables')
+        minunc = kwargs.get('minunc', 1e-10)
 
         LOGGER.debug("must finish documentation")
 
@@ -309,21 +313,22 @@ def build_matrix(self, data, sources, group=0):
         # now do a default damping target
         LOGGER.debug('damping target might be suspect for 2d objects')
         LOGGER.info("Building Damping Target")
+
         target = np.zeros(self.nunknowns, dtype=float)
         for segid, source in sources.items():
 
             for sedkey, region in source.items():
                 extid = self.sedkeys.index(sedkey)
+                if extid in self.ri:
+                    g1 = self.ri[extid]
+                    g2 = self.lamids[g1]
 
-                g1 = self.ri[extid]
-                g2 = self.lamids[g1]
-
-                if hasattr(source, 'extpars'):
-                    waves = source.extpars.wavelengths()
-                else:
-                    waves = self.defpars.wavelengths()
+                    if hasattr(source, 'extpars'):
+                        waves = source.extpars.wavelengths()
+                    else:
+                        waves = self.defpars.wavelengths()
 
-                target[g1] = region.sed(waves[g2])
+                    target[g1] = region.sed(waves[g2])
 
         self.set_damping_target(target)
 

slitlessutils/core/modules/extract/multi/multi.py

@@ -6,7 +6,7 @@
 
 from .....config import SUFFIXES, Config
 from .....logger import LOGGER
-from ....utilities import as_iterable, headers
+from ....utilities import as_iterable, get_metadata, headers
 from ...group import GroupCollection
 from ...module import Module
 from .matrix import Matrix
@@ -60,16 +60,22 @@ class Multi(Module):
     DESCRIPTION = "Extracting (multi)"
 
     def __init__(self, extorders, logdamp, mskorders=None, algorithm='golden',
-                 **kwargs):
+                 root=None, **kwargs):
         Module.__init__(self, self.extract, **kwargs, multiprocess=False)
 
+        # file root name
+        if not isinstance(root, str):
+            meta = get_metadata()
+            self.root = meta['Name']
+        else:
+            self.root = root
+
         self.extorders = as_iterable(extorders)
         self.mskorders = as_iterable(mskorders)
-
         self.optimizer = optimizer(algorithm, logdamp)
         self.matrix = Matrix(self.extorders, **kwargs)
 
-    def extract(self, data, sources, groups=None, root=None):
+    def extract(self, data, sources, groups=None):
         """
         Method to do the mult-ended spectral extraction
 
@@ -85,10 +91,6 @@ def extract(self, data, sources, groups=None, root=None):
             The collection of groups.  If set as None, then no grouping
             will be performed.  Default is None
 
-        root : str, optional
-            The root name of the output products.  If None, then
-            'slitlessutils' is used.  Default is None
-
         Notes
         -----
         This is likely not to be directly called.
@@ -116,9 +118,9 @@ def extract(self, data, sources, groups=None, root=None):
         hdul1.append(phdu)
         hdul3.append(phdu)
 
-        # file root name
-        if not isinstance(root, str):
-            root = __package__
+        # get the package meta data
+        meta = get_metadata()
+        package = meta.get('Name', '')
 
         # put loops over groups here
         if groups:
@@ -128,16 +130,16 @@ def extract(self, data, sources, groups=None, root=None):
             groups = GroupCollection()
 
         # open a PDF to write Grouping images
-        pdffile = f'{root}_{SUFFIXES["L-curve"]}.pdf'
+        pdffile = f'{self.root}_{SUFFIXES["L-curve"]}.pdf'
         LOGGER.info(f'Writing grouped L-curve figure: {pdffile}')
         with PdfPages(pdffile) as pdf:
             # add some info to the PDF
             d = pdf.infodict()
             d['Title'] = 'L-Curve Results'
             d['Author'] = getpass.getuser()
-            d['Subject'] = f'L-Curve results for grouped data from {__package__}'
-            d['Keywords'] = f'{__package__} WFSS L-curve groups'
-            d['Producer'] = __package__
+            d['Subject'] = f'L-Curve results for grouped data from {package}'
+            d['Keywords'] = f'{package} WFSS L-curve groups'
+            d['Producer'] = package
 
             for grpid, srcdict in enumerate(groups(sources)):
 
@@ -241,13 +243,13 @@ def extract(self, data, sources, groups=None, root=None):
 #
 #
 #
-#        self.matrix.lcurve.plot(f'{root}_{SUFFIXES["L-curve"]}.pdf')
+#        self.matrix.lcurve.plot(f'{self.root}_{SUFFIXES["L-curve"]}.pdf')
 
         # write out the files
         if len(hdul1) > 1:
             hdul1[0].header['FILETYPE'] = '1d spectra'
 
-            x1dfile = f'{root}_{SUFFIXES["1d spectra"]}.fits'
+            x1dfile = f'{self.root}_{SUFFIXES["1d spectra"]}.fits'
             LOGGER.info(f'Writing 1d extractions: {x1dfile}')
             hdul1.writeto(x1dfile, overwrite=True)
         else:
@@ -256,12 +258,12 @@ def extract(self, data, sources, groups=None, root=None):
         if len(hdul3) > 1:
             hdul3[0].header['FILETYPE'] = '3d spectra'
 
-            x3dfile = f'{root}_{SUFFIXES["3d spectra"]}.fits'
+            x3dfile = f'{self.root}_{SUFFIXES["3d spectra"]}.fits'
             LOGGER.info(f'Writing 3d extractions: {x3dfile}')
             hdul1.writeto(x3dfile, overwrite=True)
         else:
             LOGGER.warning('No 3d spectra written')
 
-        lcvfile = f'{root}_{SUFFIXES["L-curve"]}.fits'
+        lcvfile = f'{self.root}_{SUFFIXES["L-curve"]}.fits'
         LOGGER.info(f'Writing L-curve tabular data {lcvfile}')
         hdulL.writeto(lcvfile, overwrite=True)

slitlessutils/core/modules/extract/single/single.py

@@ -7,7 +7,7 @@
 from .....config import SUFFIXES, Config
 from .....logger import LOGGER
 from ....tables import PDTFile
-from ....utilities import headers
+from ....utilities import get_metadata, headers
 from ...module import Module
 from .boxcar import boxcar
 from .contamination import Contamination
@@ -130,7 +130,8 @@ def __init__(self, extorders, mskorders='all', savecont=False, root=None,
 
         # output file names
         if root is None:
-            root = __package__
+            root = get_metadata()['Name']
+
         self.filename = os.path.join(self.outpath, f"{root}_{SUFFIXES[self.FILETYPE]}.fits")
 
         if self.savecont:

slitlessutils/core/modules/simulate/simulate.py

@@ -1,5 +1,4 @@
 from datetime import datetime
-from importlib.metadata import metadata
 
 import numpy as np
 from astropy.io import fits
@@ -77,17 +76,13 @@ def simulate(self, data, sources, **kwargs):
            The filename of the image created
 
         """
-
         if sources.sedfile is None:
             LOGGER.critical("No SEDFile found.")
             return
 
         # grab the time
         t0 = datetime.now()
 
-        # grab the inputs
-        # insconf,insdata = data
-
         # open the table for reading
         with PDTFile(data, path=self.path, mode='r') as h5:
 
@@ -125,7 +120,6 @@ def simulate(self, data, sources, **kwargs):
 
             # process each detector
             for detname, detdata in data.items():
-
                 # load a detector
                 # detconf=detdata.config
                 # detconf=insconf[detname]
@@ -202,6 +196,7 @@ def simulate(self, data, sources, **kwargs):
                 # the function make_HDUs takes the noiseless sci, creates
                 # all ancillary data, packages into HDUs, and adds noise
                 # hdus= detdata.make_HDUs(sci,self.noisepars)
+
                 hdus = detdata.make_HDUs(sci, addnoise=self.addnoise)
 
                 # put the HDUs into the list, but first update some header info
@@ -217,18 +212,14 @@ def simulate(self, data, sources, **kwargs):
                     # append the HDU
                     hdul.append(hdu)
 
-        # record the end time
-        t1 = datetime.now()
-        dt = t1 - t0
+        # compute runtime
+        dt = datetime.now() - t0
 
         # put some times into the header
-        version = metadata(__package__).get('Version', 'unknown')
-        phdu.header.set('ORIGIN', value=f'{__package__} v{version}',
-                        after='NAXIS')
-        phdu.header.set('DATE', value=t1.strftime('%Y-%m-%d'), after='ORIGIN',
-                        comment='date this file was written (yyyy-mm-dd)')
-        phdu.header.set('RUNTIME', value=dt.total_seconds(), after='DATE',
-                        comment='run time of this file in s')
+        headers.add_preamble(phdu.header,
+                             runtime=(dt.total_seconds(), 'runtime in s'))
+        headers.add_software_log(phdu.header)
+        sources.update_header(phdu.header)    # source props
 
         # put the primary header in the HDUL
         hdul.insert(0, phdu)

slitlessutils/core/photometry/sed.py

@@ -493,9 +493,6 @@ def __call__(self, wave, fnu=False, **kwargs):
         flam = self.flam[g]
         flux = np.interp(wave, lamb, flam, **kwargs, left=flam[0], right=flam[-1])
 
-        # flux=np.interp(wave,self.lamb[g],self.flam[g],**kwargs,
-        #               left=self.flam[g[0]],right=self.flam[g[-1]])
-
         if fnu:
             flux *= ((wave / c) * (wave / 1e10))
 
@@ -626,7 +623,7 @@ def from_file(cls, filename, exten=1):
         return obj
 
     @staticmethod
-    def get_from_CDBS(atlas, filename):
+    def get_from_CDBS(atlas, filename, outpath=''):
         """
         Staticmethod to retrieve a spectrum from the Calibration Database
         System (CDBS)
@@ -639,6 +636,9 @@ def get_from_CDBS(atlas, filename):
         filename : str
             The filename in the atlas
 
+        outpath : str
+            The path for the output file.  Default is current working dir
+
         """
         # base URL for CDBS
         url = 'https://archive.stsci.edu/hlsps/reference-atlases/cdbs/grid/'
@@ -660,6 +660,13 @@ def get_from_CDBS(atlas, filename):
             return
         shutil.move(tmpfile, filename)
 
+        if isinstance(outpath, str):
+            newfile = os.path.join(outpath, filename)
+            shutil.move(filename, newfile)
+            return newfile
+        else:
+            return filename
+
     @classmethod
     def from_CDBS(obj, atlas, filename, cleanup=True):
         """

slitlessutils/core/preprocess/astrometry/list_wcs.py

@@ -29,10 +29,16 @@ def list_wcs(arg):
             obs = (h0['TELESCOP'], h0['INSTRUME'], h0['DETECTOR'])
             if obs == ('HST', 'ACS', 'WFC'):
                 filtername = h0['FILTER1']
+
+            elif obs == ('HST', 'ACS', 'SBC'):
+                filtername = h0['FILTER1']
+
             elif obs == ('HST', 'WFC3', 'UVIS'):
                 filtername = h0['FILTER']
+
             elif obs == ('HST', 'WFC3', 'IR'):
                 filtername = h0['FILTER']
+
             else:
                 raise ValueError(f"Cannot find observation {obs}")
 

slitlessutils/core/sources/source.py

@@ -117,7 +117,7 @@ def __init__(self, segid, img, seg, hdr, reg=None, zeropoint=26., grpid=0,
         self.ltv = [hdr.get('LTV1', 0.), hdr.get('LTV2', 0.)]
 
         # record some things
-        self.segid = segid  # hdr['SEGID']
+        self.segid = segid
         self.grpid = grpid
         # self.whttype = whttype
         self.backsize = max(backsize, 0)
@@ -225,6 +225,8 @@ def __init__(self, segid, img, seg, hdr, reg=None, zeropoint=26., grpid=0,
             # region image
             ri = indices.reverse(r, ignore=(0,))
             for regid, pixid in ri.items():
+                regid = int(regid)
+
                 # pixels and weights for this region
                 xx = x[pixid]
                 yy = y[pixid]
@@ -642,6 +644,8 @@ def load_sedlib(self, sedlib, throughput=None):
 
             region.sed = sed
 
+            sed.write_file(f'spectra/{self.segid}_{region.regid}.csv')
+
     def load_sed_images(self):
         pass
 
@@ -664,8 +668,6 @@ def load_sed_array(self, waves, flam):
         """
 
         for regkey, region in self.items():
-
-            # region.sed.set_sed(waves,flam*np.sum(region.w))
             region.sed.append(waves, flam * np.sum(region.w))
 
     def write_seds(self, filetype='sed', path=None, **kwargs):

slitlessutils/core/sources/sourcecollection.py

@@ -29,9 +29,9 @@ class SourceCollection(dict):
     # the default zeropoint
     DEFZERO = 26.
 
-    def __init__(self, segfile, detfile, maglim=np.inf, minpix=0, preprocessor=None,
-                 zeropoint=None, throughput=None, sedfile=None,
-                 **kwargs):
+    def __init__(self, segfile, detfile, maglim=np.inf, minpix=0,
+                 preprocessor=None, zeropoint=None, throughput=None,
+                 sedfile=None, **kwargs):
         """
         Initializer
 
@@ -196,6 +196,8 @@ def _load_classic(self, hdus, hdui, exten=0, **kwargs):
         # find pixels for each object
         ri = indices.reverse(seg, ignore=(0,))
         for segid, (y, x) in ri.items():
+            segid = int(segid)
+
             # get a bounding box
             x0 = np.maximum(np.amin(x) - self.PAD, 0)
             x1 = np.minimum(np.amax(x) + self.PAD + 1, nx - 1)

slitlessutils/core/utilities/__init__.py

@@ -2,4 +2,5 @@
 from . import headers  # noqa: F401
 from . import indices  # noqa: F401
 from .as_iterable import as_iterable  # noqa: F401
+from .get_metadata import get_metadata  # noqa: F401
 from .pool import Pool  # noqa: F401