Merge branch 'main' into develop

.readthedocs.yml.bkk

@@ -0,0 +1,15 @@
+version: 2
+
+build:
+  os: ubuntu-20.04
+  tools:
+    python: "3.10"
+
+sphinx:
+  configuration: docs/source/conf.py
+
+python:
+  install:
+    - method: pip
+      path: .
+    - requirements: docs/requirements.txt

EMToolKit/EMToolKit.py

@@ -4,12 +4,32 @@
 from ndcube import NDCube, NDCubeSequence, NDCollection
 from astropy.coordinates import SkyCoord
 from astropy.nddata import StdDevUncertainty, UnknownUncertainty
+from astropy.nddata import StdDevUncertainty, UnknownUncertainty
 from EMToolKit.schemas.basic_schemas import basic_detector, basic_source
 import astropy
 from scipy.interpolate import interp1d
 from scipy.integrate import trapezoid
 
+from scipy.interpolate import interp1d
+from scipy.integrate import trapezoid
 
+
+def em_data(maps, errs, logts, tresps, channels=None):
+    cubes = []
+    for i in range(len(maps)):
+        mapi = copy.copy(maps[i])
+        if 'CHANNEL' not in mapi.meta:
+            if channels is None:
+                mapi.meta['CHANNEL'] = mapi.meta['DETECTOR'] + mapi.meta['WAVE_STR']
+            else:
+                mapi.meta['CHANNEL'] = channels[i]
+        mapi.meta['LOGT'] = logts[i]
+        mapi.meta['TRESP'] = tresps[i]
+        mapi.meta['SHAPE'] = np.array(mapi.data.shape)
+        if 'SCHEMA' not in mapi.meta:
+            mapi.meta['SCHEMA'] = basic_detector(mapi.meta)
+        cubes.append(NDCube(mapi, uncertainty=errs[i]))
+    return NDCubeSequence(cubes)
 def em_data(maps, errs, logts, tresps, channels=None):
     cubes = []
     for i in range(len(maps)):
@@ -39,6 +59,27 @@ def dem_model(coeffs, logts, bases, coord_info, algorithm, wrapper, meta=None, w
             print('Warning in EMToolKit dem_model: coord_info is not a wcs or dict')
         wcs = meta.get('wcs', None)
 
+    if meta is None:
+        if wcs is not None:
+            meta = dict(wcs.to_header())
+        else:
+            print('Warning in EMToolKit dem_model: need wcs or image meta')
+    if 'LOGT' not in meta:
+        meta['LOGT'] = logts[0]
+    if 'SCHEMA' not in meta:
+        meta['SCHEMA'] = basic_source([Map(coeffs[0], meta)])
+    for i in range(nd):
+        logt0 = np.median(logts[i][np.where(bases[i] == np.max(bases[i]))])
+        metai = {
+            'BASIS': bases[i],
+            'LOGT': logts[i],
+            'LOGT0': logt0,
+            'ALGORITHM': algorithm,
+            'WRAPPER': wrapper,
+            'WRAPARGS': wrapargs,
+        }
+        dem_sequence.append(NDCube(coeffs[i], wcs=wcs, meta={**meta, **metai}))
+    return NDCubeSequence(dem_sequence, meta={'ALGORITHM': algorithm, 'WRAPPER': wrapper, 'WRAPARGS': wrapargs})
     if meta is None:
         if wcs is not None:
             meta = dict(wcs.to_header())
@@ -68,6 +109,12 @@ def __init__(self, datasequence):
         self.collection = NDCollection([("data", datasequence), ("models", [])])
         self.precomputed_interpolations = None
 
+    def data(self):
+        return self.collection['data']
+    def __init__(self, datasequence):
+        self.collection = NDCollection([("data", datasequence), ("models", [])])
+        self.precomputed_interpolations = None
+
     def data(self):
         return self.collection['data']
 
@@ -147,6 +194,33 @@ def synthesize_data(self, logts, tresps, algorithm=None, channels=None, ilo=0, i
         [synthmaps, syntherrs] = [[], []]
         self.precompute_interpolations()
 
+        for i in range(ndata):
+            synthdata = np.zeros(model[0].data[ilo:ihi, jlo:jhi].shape)
+            syntherrs.append(UnknownUncertainty(np.zeros(model[0].data.shape) - 1))
+            synthmap = copy.deepcopy(model[0])[ilo:ihi, jlo:jhi]
+            synthmap.meta['NAXIS'] = 2
+            synthmap.meta['ALGORITHM'] = algorithm
+            synthmap.meta['CHANNEL'] = channels[i]
+            datainterp = interp1d(logts[i], tresps[i], fill_value=0.0, bounds_error=False)
+            for ind, component in enumerate(model):
+                basisinterp = self.precomputed_interpolations[algorithm][ind] #[model.index(component)]
+                logt = np.unique(np.hstack([component.meta['LOGT'], logts[i]]))
+                coupling = trapezoid(datainterp(logt) * basisinterp(logt), x=logt)
+                synthdata += coupling * component.data[ilo:ihi, jlo:jhi]
+            synthmap.data[:] = synthdata[:]
+            synthmaps.append(synthmap)
+    def synthesize_data(self, logts, tresps, algorithm=None, channels=None, ilo=0, ihi=-1, jlo=0, jhi=-1, meta=None):
+        if logts[0].size == 1:
+            [logts, tresps] = [[logts], [tresps]]
+        ndata = len(logts)
+        if algorithm is None:
+            algorithm = self.collection['models'][0]
+        if channels is None:
+            channels = ['SYNTHDATA' + str(i) for i in range(ndata)]
+        model = self.collection[algorithm]
+        [synthmaps, syntherrs] = [[], []]
+        self.precompute_interpolations()
+
         for i in range(ndata):
             synthdata = np.zeros(model[0].data[ilo:ihi, jlo:jhi].shape)
             syntherrs.append(UnknownUncertainty(np.zeros(model[0].data.shape) - 1))
@@ -164,6 +238,7 @@ def synthesize_data(self, logts, tresps, algorithm=None, channels=None, ilo=0, i
             synthmaps.append(synthmap)
 
         return em_data(synthmaps, syntherrs, logts, tresps, channels=channels)
+        return em_data(synthmaps, syntherrs, logts, tresps, channels=channels)
 
     def synthesize_map(self, map, logt=None, tresp=None, algorithm=None, channel=None):
         if 'CHANNEL' not in map.meta:

EMToolKit/algorithms/simple_reg_dem_wrapper.py

@@ -10,6 +10,36 @@
 import EMToolKit.EMToolKit as emtk
 
 
+def simple_reg_dem_wrapper(datasequence, wrapargs=None):
+    """
+    Wrapper function for the simple regularized Differential Emission Measure (DEM) calculation.
+
+    This function prepares input data and passes it to the `simple_reg_dem` algorithm.
+    It processes the input data to ensure that all input maps have consistent dimensions,
+    extracts the necessary metadata, and then calls the DEM calculation.
+
+    Parameters
+    ----------
+    datasequence : NDCubeSequence
+        A sequence of data cubes containing the observations. Each cube should contain
+        2D spatial data with associated uncertainties and metadata.
+    wrapargs : dict, optional
+        Additional arguments to pass to the initialization routines of the `simple_reg_dem` function.
+
+    Returns
+    -------
+    list of numpy.ndarray
+        The calculated DEM coefficients for each temperature bin.
+    list of numpy.ndarray
+        The temperature bins used in the calculation.
+    list of numpy.ndarray
+        The basis functions for the temperature bins.
+    WCS
+        World Coordinate System (WCS) information from the input data.
+    str
+        A string identifier for the DEM method used.
+    """
+    # Determine the smallest data cube dimensions in the sequence
 def simple_reg_dem_wrapper(datasequence, wrapargs=None):
     """
     Wrapper function for the simple regularized Differential Emission Measure (DEM) calculation.
@@ -51,6 +81,21 @@ def simple_reg_dem_wrapper(datasequence, wrapargs=None):
     tresps = np.zeros([logt.size, nc])  # Temperature response functions
     exptimes = np.zeros(nc)  # Exposure times
 
+    # Fill the data, uncertainty, and metadata arrays
+    for i in range(nc):
+        datacube[:, :, i] = datasequence[i].data[0:nx, 0:ny]
+        errscube[:, :, i] = datasequence[i].uncertainty.array[0:nx, 0:ny]
+        tresps[:, i] = datasequence[i].meta['tresp']
+    [nx, ny] = datasequence[0].data.shape
+    for seq in datasequence:
+        [nx, ny] = [np.min([seq.data.shape[0], nx]), np.min([seq.data.shape[1], ny])]
+
+    logt = datasequence[0].meta['logt']  # Temperature bins
+    datacube = np.zeros([nx, ny, nc])  # Data cube for storing observations
+    errscube = np.zeros([nx, ny, nc])  # Data cube for storing uncertainties
+    tresps = np.zeros([logt.size, nc])  # Temperature response functions
+    exptimes = np.zeros(nc)  # Exposure times
+
     # Fill the data, uncertainty, and metadata arrays
     for i in range(nc):
         datacube[:, :, i] = datasequence[i].data[0:nx, 0:ny]
@@ -61,17 +106,31 @@ def simple_reg_dem_wrapper(datasequence, wrapargs=None):
     # Perform the DEM calculation
     coeffs, chi2 = simple_reg_dem(datacube, errscube, exptimes, logt, tresps)
 
+    # Perform the DEM calculation
+    coeffs, chi2 = simple_reg_dem(datacube, errscube, exptimes, logt, tresps)
+
     # Simple_reg_dem puts the temperature axis last. Transpose so it's the first:
     coeffs = coeffs.transpose(np.roll(np.arange(coeffs.ndim), 1))
+    coeffs = coeffs.transpose(np.roll(np.arange(coeffs.ndim), 1))
 
     nt = logt.size
     wcs = datasequence[0].wcs  # WCS information from the first cube
     basislogt = np.linspace(np.min(logt), np.max(logt), 2 * (nt - 1) + 1)
     logts, bases = [], []
 
+    # Create the basis functions for the temperature bins
+    for i in range(nt):
+    wcs = datasequence[0].wcs  # WCS information from the first cube
+    basislogt = np.linspace(np.min(logt), np.max(logt), 2 * (nt - 1) + 1)
+    logts, bases = [], []
+
     # Create the basis functions for the temperature bins
     for i in range(nt):
         basis = (basislogt == logt[i]).astype(np.float64)
+        if i > 0:
+            basis += (basislogt - logt[i - 1]) * (basislogt < logt[i]) * (basislogt > logt[i - 1]) / (logt[i] - logt[i - 1])
+        if i < nt - 1:
+            basis += (logt[i + 1] - basislogt) * (basislogt < logt[i + 1]) * (basislogt > logt[i]) / (logt[i + 1] - logt[i])
         if i > 0:
             basis += (basislogt - logt[i - 1]) * (basislogt < logt[i]) * (basislogt > logt[i - 1]) / (logt[i] - logt[i - 1])
         if i < nt - 1:
@@ -124,10 +183,14 @@ def autoloading_simple_reg_dem_wrapper(datasequence, data_dir=".data/default/",
         tstart = time.time()
         simpl_out = simple_reg_dem_wrapper(datasequence, wrapargs)
         print('Done! Simple method took', time.time() - tstart)
+        print('Done! Simple method took', time.time() - tstart)
         simple_reg_demsequence = emtk.dem_model(*simpl_out, simple_reg_dem_wrapper)
 
+        # Save the DEM sequence to a file
+
         # Save the DEM sequence to a file
         with open(pk_file, 'wb') as file:
             pickle.dump((simple_reg_demsequence, simpl_out), file)
 
+
     return simple_reg_demsequence, simpl_out

EMToolKit/algorithms/sparse_em_wrapper.py

@@ -64,6 +64,12 @@ def sparse_em_wrapper(datasequence, wrapargs={}):
     Dict, lgtaxis, basis_funcs, bases_sigmas = sparse_em_init(trlogts, tresps, differential=True)
     coeffs, zmax, status = sparse_em_solve(datacube, errscube, exptimes, Dict)
 
+    # Sparse_em_solve puts the temperature axis last. Transpose so it's the first:
+    coeffs = coeffs.transpose(np.roll(np.arange(coeffs.ndim), 1))
+    # Initialize and solve the sparse EM problem
+    Dict, lgtaxis, basis_funcs, bases_sigmas = sparse_em_init(trlogts, tresps, differential=True)
+    coeffs, zmax, status = sparse_em_solve(datacube, errscube, exptimes, Dict)
+
     # Sparse_em_solve puts the temperature axis last. Transpose so it's the first:
     coeffs = coeffs.transpose(np.roll(np.arange(coeffs.ndim), 1))
 

EMToolKit/instruments/xrt.py

@@ -78,7 +78,6 @@ def xrt_temperature_response(map_in, temperature_array, *, do_plot=False):
 
 	channel = map_in.meta['ec_fw1_']+'/'+map_in.meta['ec_fw2_']
 	channel = channel.replace("Open/", "")
-
 	filt = channel
 	date_time = "2007-09-22T21:59:59"
 	Temperature_Response_Fundamental = xrtpy.response.TemperatureResponseFundamental(

docs/requirements.txt

@@ -2,6 +2,32 @@ sphinx
 sphinx-rtd-theme
 sphinx-autobuild
 sphinx-autodoc-typehints
+sphinx-gallery
+sphinx-hoverxref
+sphinx-toggleprompt
+sphinx-copybutton
+sphinx-design
+nbsphinx
+nbconvert
+nbsphinx-link
+ipython
+ipywidgets
+pandoc==2.3
+matplotlib
+myst-parser
+esbonio>=0.12.0
+
+sunpy[all]
+astropy
+scipy
+numba
+tqdm
+ndcube
+lxml
+zeep
+drms
+memory-profiler
+
 sphinx-automodapi
 sphinx-gallery
 sphinx-hoverxref

docs/source/EMToolKit.algorithms.rst

@@ -0,0 +1,82 @@
+EMToolKit.algorithms package
+============================
+
+EMToolKit.algorithms.multi\_inst\_dem\_wrapper module
+-----------------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.multi_inst_dem_wrapper
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.simple\_reg\_dem module
+--------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.simple_reg_dem
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.simple\_reg\_dem\_parallel module
+------------------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.simple_reg_dem_parallel
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.simple\_reg\_dem\_wrapper module
+-----------------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.simple_reg_dem_wrapper
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.simple\_reg\_dem\_wrapper\_parallel module
+---------------------------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.simple_reg_dem_wrapper_parallel
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.sparse\_em module
+--------------------------------------
+
+.. automodule:: EMToolKit.algorithms.sparse_em
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.sparse\_em\_wrapper module
+-----------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.sparse_em_wrapper
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.sparse\_nlmap\_dem\_wrapper module
+-------------------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.sparse_nlmap_dem_wrapper
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.algorithms.sparse\_nlmap\_solver module
+-------------------------------------------------
+
+.. automodule:: EMToolKit.algorithms.sparse_nlmap_solver
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Module contents
+---------------
+
+.. automodule:: EMToolKit.algorithms
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/source/EMToolKit.instruments.rst

@@ -0,0 +1,18 @@
+EMToolKit.instruments package
+=============================
+
+EMToolKit.instruments.aia module
+--------------------------------
+
+.. automodule:: EMToolKit.instruments.aia
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+EMToolKit.instruments.xrt module
+--------------------------------
+
+.. automodule:: EMToolKit.instruments.xrt
+   :members:
+   :undoc-members:
+   :show-inheritance: