adds metadata to output

.gitignore

@@ -165,3 +165,4 @@ data/
 output/
 *.fits
 */*.fits
+!tests/*.fits

.pre-commit-config.yaml

@@ -26,6 +26,7 @@ repos:
       - id: check-yaml
       - id: debug-statements
       - id: check-added-large-files
+        args: [ '--enforce-all','--maxkb=1000' ]
       - id: end-of-file-fixer
         exclude: ".*(.fits|.fts|.fit|.txt|.pro|.asdf|.bib|tca.*)"
       - id: mixed-line-ending

docs/configuration.rst

@@ -23,7 +23,7 @@ is divided into sections:
 There are four configurables for this section:
 
 - *overlappogram*: path to the overlappogram image to be inverted.
-- *weights*: path to the accompanying weights used in the inversion.
+- *weights*: path to the accompanying weights used in the inversion. Weights should be in units of :math:`\frac{1}{\sigma}` where :math:`\sigma` is the uncertainty or standard deviation
 - *response*: path to the instrument response.
 - *gnt*: path to the file containing atomic physics values from Chianti, the *G(n, t)* function.
 

example_config.toml

@@ -19,7 +19,7 @@ response_dependency_list = [5.7, 5.8, 5.9, 6.0 , 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6
 smooth_over = 'dependence'
 
 [model]
-alphas = [3E-5, 4E-5, 0.1] #[0.2, 0.1, 0.01, 0.005]
+alphas = [3E-5] #, 4E-5, 0.1] #[0.2, 0.1, 0.01, 0.005]
 rhos = [0.1]
 warm_start = false
 tol = 1E-2

overlappogram/cli.py

@@ -50,13 +50,13 @@ def unfold(config):
             print(f"Beginning inversion for alpha={alpha}, rho={rho}.")
             start = time.time()
             em_cube, prediction, scores, unconverged_rows = inversion.invert(
-                overlappogram,
-                config["model"],
-                alpha,
-                rho,
-                num_threads=config["execution"]["num_threads"],
-                mode_switch_thread_count=config["execution"]["mode_switch_thread_count"],
-                mode=MODE_MAPPING.get(config['execution']['mode'], 'invalid')
+                    overlappogram,
+                    config["model"],
+                    alpha,
+                    rho,
+                    num_threads=config["execution"]["num_threads"],
+                    mode_switch_thread_count=config["execution"]["mode_switch_thread_count"],
+                    mode=MODE_MAPPING.get(config['execution']['mode'], 'invalid')
             )
             end = time.time()
             print(
@@ -70,13 +70,15 @@ def unfold(config):
             )
             save_em_cube(
                 em_cube,
-                os.path.join(config["output"]["directory"], f"{config['output']['prefix']}_emcube_{postfix}.fits"),
+                os.path.join(config["output"]["directory"],
+                             f"{config['output']['prefix']}_emcube_{postfix}.fits"),
                 config["output"]["overwrite"],
             )
 
             save_prediction(
                 prediction,
-                os.path.join(config["output"]["directory"], f"{config['output']['prefix']}_prediction_{postfix}.fits"),
+                os.path.join(config["output"]["directory"],
+                             f"{config['output']['prefix']}_prediction_{postfix}.fits"),
                 overwrite=config["output"]["overwrite"],
             )
 

overlappogram/error.py

@@ -4,3 +4,15 @@ class OverlappogramWarning(Warning):
 
 class NoWeightsWarnings(OverlappogramWarning):
     """There are no weights passed to unfold."""
+
+
+class OverlappogramError(Exception):
+    pass
+
+
+class InvalidDataFormatError(OverlappogramError):
+    """The data file has an error that prevents processing."""
+
+
+class InvalidInversionModeError(OverlappogramError):
+    """The inversion mode is not allowed."""

overlappogram/inversion.py

@@ -5,13 +5,14 @@
 from enum import Enum
 from threading import Lock
 
+import astropy.wcs as wcs
 import numpy as np
 from ndcube import NDCube
 from sklearn.exceptions import ConvergenceWarning
 from sklearn.linear_model import ElasticNet
 from tqdm import tqdm
 
-from overlappogram.error import NoWeightsWarnings
+from overlappogram.error import InvalidInversionModeError, NoWeightsWarnings
 from overlappogram.response import prepare_response_function
 
 __all__ = ["Inverter"]
@@ -47,8 +48,8 @@ def __init__(
         self._unconverged_rows = []
 
         self._overlappogram: NDCube | None = None
-        self._em_data: NDCube | None = None
-        self._inversion_prediction: NDCube | None = None
+        self._em_data: np.ndarray | None = None
+        self._inversion_prediction: np.ndarray | None = None
         self._row_scores: np.ndarray | None = None
         self._overlappogram_width: int | None = None
         self._overlappogram_height: int | None = None
@@ -61,6 +62,11 @@ def __init__(
             field_angle_range=field_angle_range,
             response_dependency_list=response_dependency_list,
         )
+        self._response_meta = response_cube.meta
+
+        if response_dependency_list is not None:
+            self._response_meta['temperatures'] = np.array([(i, t) for i, t in enumerate(response_dependency_list)],
+                                                           dtype=[('index', '>i2'), ('logt', '>f4')])
 
         self._progress_bar = None  # initialized in invert call
 
@@ -111,7 +117,6 @@ def _progress_indicator(self, future):
         with self._thread_count_lock:
             if not future.cancelled():
                 self._completed_row_count += 1
-                #print(f"{self._completed_row_count / self.total_row_count * 100:3.0f}% complete", end="\r")
                 self._progress_bar.update(1)
 
     def _switch_to_row_inversion(self, model_config, alpha, rho, num_row_threads=50):
@@ -249,29 +254,6 @@ def _initialize_with_overlappogram(self, overlappogram):
         self._inversion_prediction = np.zeros((self._overlappogram_height, self._overlappogram_width), dtype=np.float32)
         self._row_scores = np.zeros((self._overlappogram_height, 1), dtype=np.float32)
 
-        # TODO : add metadata to outputs
-        # self.inv_date = (
-        #     datetime.datetime.now()
-        #     .isoformat(timespec="milliseconds")
-        #     .replace("+00:00", "Z")
-        # )
-        # # try:
-        #     self.rsp_func_date = rsp_func_hdul[0].header["DATE"]
-        # except KeyError:
-        #     self.rsp_func_date = ""
-        # try:
-        #     self.abundance = rsp_func_hdul[0].header["ABUNDANC"]
-        # except KeyError:
-        #     self.abundance = ""
-        # try:
-        #     self.electron_distribution = rsp_func_hdul[0].header["ELECDIST"]
-        # except KeyError:
-        #     self.electron_distribution = ""
-        # try:
-        #     self.chianti_version = rsp_func_hdul[0].header["CHIANT_V"]
-        # except KeyError:
-        #     self.chianti_version = ""
-
     def invert(
         self,
         overlappogram: NDCube,
@@ -305,16 +287,18 @@ def invert(
             self._start_row_inversion(model_config, alpha, rho, num_threads)
             self._collect_results(np.inf, model_config, alpha, rho)
         else:
-            raise ValueError("Invalid InversionMode.")
+            raise InvalidInversionModeError("Invalid InversionMode.")
 
         for executor in self.executors:
             executor.shutdown()
 
         self._progress_bar.close()
 
+        out_wcs = wcs.WCS(naxis=2)
+
         return (
-            np.transpose(self._em_data, (2, 0, 1)),
-            self._inversion_prediction,
+            NDCube(data=np.transpose(self._em_data, (2, 0, 1)), wcs=out_wcs, meta=self._response_meta),
+            NDCube(data=self._inversion_prediction, wcs=out_wcs, meta=self._response_meta),
             self._row_scores,
             self._unconverged_rows,
         )

overlappogram/io.py

@@ -1,11 +1,24 @@
+from datetime import datetime
+
 from astropy.io import fits
 from astropy.nddata import StdDevUncertainty
 from astropy.wcs import WCS
 from ndcube import NDCube
 
+from .error import InvalidDataFormatError
+
 __all__ = ["load_overlappogram", "load_response_cube", "save_em_cube", "save_spectral_cube", "save_prediction"]
 
 
+RESPONSE_HEADER_KEYS = ['DATE',
+                        'VALUE',
+                        'FIELDANG',
+                        'RSP_DATE',
+                        'DEPNAME',
+                        'ABUNDANC',
+                        'ELECDIST']
+
+
 def load_overlappogram(image_path: str, weights_path: str | None) -> NDCube:
     with fits.open(image_path) as image_hdul:
         image = image_hdul[0].data
@@ -31,16 +44,52 @@ def load_response_cube(path: str) -> NDCube:
         field_angles = hdul[2].data
     meta = dict(header)
     meta.update({"temperatures": temperatures, "field_angles": field_angles})
+    for key in RESPONSE_HEADER_KEYS:
+        meta[key] = header[key]
     return NDCube(response, wcs=wcs, meta=meta)
 
 
-def save_em_cube(cube, path: str, overwrite: bool = True) -> None:
-    fits.writeto(path, cube, overwrite=overwrite)
+def save_em_cube(cube: NDCube, path: str, overwrite: bool = True) -> None:
+    if "temperatures" not in cube.meta:
+        raise InvalidDataFormatError("Temperatures are missing form the EM cube's metadata.")
+
+    hdul = fits.HDUList([fits.PrimaryHDU(cube.data),
+                         fits.BinTableHDU(cube.meta['temperatures']),
+                         # todo : make sure this is the right temperatures and not just a pass through
+                         ])
+    for key in RESPONSE_HEADER_KEYS:
+        hdul[0].header[key] = cube.meta[key]
+    hdul[0].header['DATE'] = datetime.now().isoformat()
+    hdul.writeto(path, overwrite=overwrite)
+
+
+def save_prediction(prediction: NDCube, path: str, overwrite: bool = True) -> None:
+    if "temperatures" not in prediction.meta:
+        raise InvalidDataFormatError("Temperatures are missing form the prediction cube's metadata.")
+
+    hdul = fits.HDUList([fits.PrimaryHDU(prediction.data),
+                         fits.BinTableHDU(prediction.meta['temperatures']),
+                         # todo : make sure this is the right temperatures and not just a pass through
+                         ])
+    for key in RESPONSE_HEADER_KEYS:
+        hdul[0].header[key] = prediction.meta[key]
+    hdul[0].header['DATE'] = datetime.now().isoformat()
+    hdul.writeto(path, overwrite=overwrite)
 
 
-def save_prediction(prediction, path: str, overwrite: bool = True) -> None:
-    fits.writeto(path, prediction, overwrite=overwrite)
+def save_spectral_cube(spectral_cube: NDCube, path: str, overwrite: bool = True) -> None:
+    if "temperatures" not in spectral_cube.meta:
+        raise InvalidDataFormatError("Temperatures are missing form the spectral cube's metadata.")
 
+    if "ions" not in spectral_cube.meta:
+        raise InvalidDataFormatError("Ions are missing form the spectral cube's metadata.")
 
-def save_spectral_cube(spectral_cube, path: str, overwrite: bool = True) -> None:
-    fits.writeto(path, spectral_cube, overwrite=overwrite)
+    hdul = fits.HDUList([fits.PrimaryHDU(spectral_cube.data),
+                         fits.BinTableHDU(spectral_cube.meta['temperatures']),
+                         fits.BinTableHDU(spectral_cube.meta['ions']),
+                         # todo : make sure this is the right temperatures and not just a pass through
+                         ])
+    for key in RESPONSE_HEADER_KEYS:
+        hdul[0].header[key] = spectral_cube.meta[key]
+    hdul[0].header['DATE'] = datetime.now().isoformat()
+    hdul.writeto(path, overwrite=overwrite)

overlappogram/response.py

@@ -7,6 +7,8 @@
 def prepare_response_function(
     response_cube: NDCube, field_angle_range=None, response_dependency_list=None, fov_width=2
 ) -> (np.ndarray, float, float):
+    # from Dyana Beabout
+
     num_dep, num_field_angles, rsp_func_width = np.shape(response_cube.data)
 
     dependency_list = [t for (_, t) in response_cube.meta["temperatures"]]

overlappogram/spectral.py

@@ -1,17 +1,22 @@
 
+import astropy.wcs as wcs
 import numpy as np
 from astropy.io import fits
+from ndcube import NDCube
 
 __all__ = [
     "create_spectrally_pure_images",
 ]
 
 
-def create_spectrally_pure_images(image_list: list, gnt_path: str, rsp_dep_list: list):
-    # Create output directory.
+def create_spectrally_pure_images(image_list: list[NDCube],
+                                  gnt_path: str,
+                                  rsp_dep_list: list | None) -> NDCube:
+    # from Dyana Beabout
     num_images = len(image_list)
     if num_images > 0:
         with fits.open(gnt_path) as gnt_hdul:
+            ions = gnt_hdul[2].data
             gnt_data_values = gnt_hdul[0].data.astype(np.float64)
             num_gnts, num_gnt_deps = np.shape(gnt_data_values)
             gnt_dep_list = gnt_hdul[1].data["logt"]
@@ -47,7 +52,7 @@ def create_spectrally_pure_images(image_list: list, gnt_path: str, rsp_dep_list:
             for index in range(len(image_list)):
                 # Create spectrally pure data cube.
                 for em_data in image_list:
-                    em_data_cube = em_data.astype(np.float64)
+                    em_data_cube = em_data.data.astype(np.float64)
                     em_data_cube = np.transpose(em_data_cube, axes=(1, 2, 0))
                     if index == 0:
                         image_height, num_slits, num_logts = np.shape(em_data_cube)
@@ -60,4 +65,13 @@ def create_spectrally_pure_images(image_list: list, gnt_path: str, rsp_dep_list:
                             em_data_cube[:, :, 0:num_rsp_deps] * 10**26 * gnt_values[gnt_num, 0:num_rsp_deps]
                         ).sum(axis=2)
                         gnt_data_cube[:, :, gnt_num] = gnt_image
-    return np.transpose(gnt_data_cube, (2, 0, 1))
+
+    out_wcs = wcs.WCS(naxis=2)
+    out = NDCube(data=np.transpose(gnt_data_cube, (2, 0, 1)),
+                 wcs=out_wcs,
+                 meta={
+                     "temperatures": image_list[0].meta["temperatures"],
+                     "ions": ions})
+    out.meta.update(image_list[0].meta)
+
+    return out

tests/test_config.toml

@@ -0,0 +1,32 @@
+[paths]
+overlappogram = "./tests/test_overlappogram.fits"
+weights = "./tests/test_weights.fits"
+response = "./tests/test_response.fits"
+gnt = "data/ECCCO_speedtest_runs/master_gnt_eccco_inphotons_cm3persperpix_with_tables.fits"
+
+[output]
+prefix = "test"
+make_spectral = true
+directory = "output/example/"
+overwrite = true
+
+[inversion]
+solution_fov_width = 2
+detector_row_range = [0, 9]
+field_angle_range = [-1227, 1227]
+response_dependency_name = "logt"
+response_dependency_list = [5.7, 5.8, 5.9, 6.0 , 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8]
+smooth_over = 'dependence'
+
+[model]
+alphas = [3E-5] #, 4E-5, 0.1] #[0.2, 0.1, 0.01, 0.005]
+rhos = [0.1]
+warm_start = false
+tol = 1E-2
+max_iter = 1_000
+selection = 'cyclic'
+
+[execution]
+num_threads = 2
+mode_switch_thread_count = 100  # only used in hybrid mode
+mode = "row"