Implement regularized inversion method of Hannah and Kontar (2012)

[wtbarnes]

The regularized inversion method of Hannah and Kontar (2012) is one of the most popular and widely used DEM inversion methods, particularly with AIA data.

There are both Python and IDL implementations:

• Python
• Alternative Python
• IDL

From discussions with @alasdairwilson, it seems that the Python version is extremely fast and matches the output of the original IDL code quite well and so would be an excellent candidate for trying out in this framework.

[PaulJWright]

I have a write-up of the code that could be a nice starting point as documentation for this.

[wtbarnes]

That'd be great. We should think about how we want to structure the documentation and this would be a nice test case for that

[PaulJWright]

Yes, if we can figure out a docs structure I can start looking at adding things.

[wtbarnes]

I created #20 to have this discussion since the scope is outside of just this method.

[wtbarnes]

I used @alasdairwilson's version of demregpy and actually got a version working:

class HK12Model(GenericModel):
    
    def _model(self, alpha=1.0, increase_alpha=1.5, max_iterations=10, guess=None, return_em=False):
        errors = np.array([c.uncertainty.array.squeeze() for c in self.data]).T
        dem,edem,elogt,chisq,dn_reg=dn2dem_pos(
            self.data_matrix.value.T,
            errors,
            self.kernel_matrix.value.T,
            np.log10(self.temperature_bin_centers.to(u.K).value),
            self.temperature_bin_edges.to(u.K).value,
            max_iter=max_iterations,
            reg_tweak=alpha,
            rgt_fact=increase_alpha,
            dem_norm0=guess,
        )
        dem_unit = self.data_matrix.unit / self.kernel_matrix.unit / self.temperature_bin_edges.unit
        dem = dem * dem_unit
        edem = edem * dem_unit
        if return_em:
            delta_t = np.diff(self.temperature_bin_edges)
            dem = dem * delta_t
            edem = edem * delta_t
        return dem.T, edem.T
        
    @classmethod
    def defines_model_for(self, *args, **kwargs):
        return kwargs.get('model') == 'hk12'

Ideally this would of course be a PR, but there are some things we should think about:

• Should demregpy just be implemented directly in this package? Or should it be a separate package itself? At the moment, it is difficult to depend on this code externally as it is not packaged and available via pip or conda.
• This code also returns errors in temperature space. Currently, sunkit-dem only has the ability to return errors in the DEM
• It would be nice to also return the chi-squared values as well, though again, there's not really functionality for this yet. Maybe we should be returning an NDCollection rather than an NDCube so we can return the temperature errors and the chi-squared maps which have the same axes in space.

[PaulJWright]

• Should demregpy just be implemented directly in this package? Or should it be a separate package itself? At the moment, it is difficult to depend on this code externally as it is not packaged and available via pip or conda.

I think it makes sense to have demregpy packaged itself, and a dependency of this package. That however does then shift responsibility to @alasdairwilson or @ianan to maintain it; on the otherhand, having it implemented directly in this package also provides a good template of how to implement other methods.

• This code also returns errors in temperature space. Currently, sunkit-dem only has the ability to return errors in the DEM

Yes, I think it should have that functionality

• It would be nice to also return the chi-squared values as well, though again, there's not really functionality for this yet. Maybe we should be returning an NDCollection rather than an NDCube so we can return the temperature errors and the chi-squared maps which have the same axes in space.

Yes, agreed.

[PaulJWright]

I can make this into its own issue (unsure if it's a sunkit-dem or NDCube issue), but just an FYI for now.

When implementing this method with NDCube 2.0.0 I run into issues with NDSequence for 1D and 0D data. Specifically, the issue is here, and is a direct result of the following NDCube behaviour.

sunkit-dem/sunkit_dem/base_model.py

Line 153 in 9a439c9

wcs = self.data[0].wcs.to_header()

For the following 3D sequence seq, the behaviour is as expected:

>>> seq.data[0].dimensions

[1, 50, 50] pix

>>> seq.data[0].wcs

WCSAXES =                    3 / Number of coordinate axes                      
CRPIX1  =      35.372013038278 / Pixel coordinate of reference point            
CRPIX2  =       65.40190861244 / Pixel coordinate of reference point            
CRPIX3  =                  1.0 / Pixel coordinate of reference point     
...

where the type(seq.data[0].wcs) is astropy.wcs.wcs.WCS.

In the 1D and 0D cases, the following is observed:

>>> slice_1d = seq.index_as_cube[0:6,pix_lat,:]
>>> slice_1d.data[0].dimensions
[1, 50] pix

>>> slice_1d.data[0].wcs.to_header()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-55-863c7d91680f> in <module>
----> 1 slice_0d.data[0].wcs.to_header()

AttributeError: 'HighLevelWCSWrapper' object has no attribute 'to_header'

as type(slice_1d.data[0].wcs) is astropy.wcs.wcsapi.high_level_wcs_wrapper.HighLevelWCSWrapper

Unsure if there is an easy workaround here, but thought I'd mention before opening an issue

[wtbarnes]

Thanks for pointing this out. This is constantly an issue when dealing with sliced WCSs. Now that ndcube 2.0 is out and stable, we should really revisit all of that code (that I originally wrote against ndcube 1.x).

[PaulJWright]

Thanks for pointing this out. This is constantly an issue when dealing with sliced WCSs. Now that ndcube 2.0 is out and stable, we should really revisit all of that code (that I originally wrote against ndcube 1.x).

Are we now in a position to integrate this method?