Refactor of Fitting Infrastructure

[DanRyanIrish]

Provide a general description of the issue or problem.

The fitting infrastructure provided by @KriSun95 is a valuable addition to the package. To help modularity, flexibility and ease of maintenance, a refactor would be beneficial. This issue sketches out how this could be done.

Components

• Data structure: A structure of holding the spectrum/spectra. (Inherits from NDCube?/xarray?/Spectrum1D?)
• Photon models: Functions providing photons as a function of energy determined by certain physical variables, e.g. temperature (Thermal and thin/thick target models exist but need more work.)
• Spectral Response Model: Given physical variables of photon model/composition of photon models, returns count spectrum. It is a pipeline of composed of:
  • photon model/photons models
  • instrument response including but not limited to detector response matrix
• Fitting tool: Compares observed and model counts and iterates physical variables to minimise residuals.
• These tools must be able to handle multiple spectra from multiple instruments with different spectral responses and minimise fit using the same photon model(s).

Data Structure(s) Requirements

• Stores and provides easy access to count data and spectral axis values for a single spectrum from a single instrument.
• May also require to hold multiple spectra from a single instrument.
• A collection object must store and count and spectral axis values for spectra from multiple instruments.
• Handles irregular energy bin widths
• Can return energy bin centers, edges and widths.

Spectral Response Model

• Photon model/composition of photon models
• Instrumental response model
  • Can be a fixed matrix or have instrument variables, e.g. thickness of solar black
• Free and fixed parameters to photon models and possible instrumental variables

Spectral Response Model Collection

• Take one photon model/composition photon model and pass them to multiple Spectral Response Models, e.g. different instruments.

• More text to follow...

[settwi]

I like the general setup of this. I think we should avoid inheritance unless it will offer an obvious benefit to the code. We can just use composition (i.e. a "has-a" NDCube rather than an "is-a" relationship). It makes it easier to swap things out later.

@KriSun95 have been thinking of having a mini-workshop for a few days at some point in the near future to hash out some ideas too.

[samaloney]

I think the main thing is to define the API. How that API is implemented e.g composition, inheritance or a mix is a detail, all be it an important one. Ideally we would have something something like the NDCube SEP realistically probably a simlified version 😆

[DanRyanIrish]

After the discussion at the STIX meeting last week, here's an update on one option for the refactor

Overall structure for fitting a single model to a spectrum/multiple spectra on same spectral grid

3 main objects:

• Spectrum: Holds the observed data.
  • Suggest this follows the API of NDCube and any extra methods/attributes follows the API of specutils.Spectrum1D (which is an instance of NDCube).
  • A spectrum object should be N-D, so as to store multiple spectra on the same spectral grid. This should facilitate simultaneous vectorised fitting of multiple spectra with the same model.
  • The reason for not using Spectrum1D directly is that it provides a lot of additional attributes and methods that are not needed in solar X-ray spectral fitting
  • NDCube.mask can be used to exclude certain spectral areas from the fitting.
• Model: The model to fit to the spectrum/spectra
  • Suggest using the astropy.modeling API for this, even if certain model classes must be implemented from scratch.
  • This prevents us defining yet another API, will make it easier for our infrastructure to be compatible with the astropy.modeling framework, and be more intuitive to users already used to this framework.
  • This will allow us to use astropy's CompoundModel API where each stage in the model chain is represented by a submodel, e.g. model = (thermal_model + thick_target_model) | optical_path_model | drm_model
• Fitter: Fits the data to the model.
  • Could use astropy.fitting or could be a bespoke fitter. More likely the latter.

What's Needed

Fitter

• Confirm that the astropy.modeling API enables vectorised parameters and whether there is any problem with using the API in our own fitting framework

Model

• All model components, (e.g. photon model, effective area, attenuator, etc.) should be possible to describe with astropy models, or at least something bespoke that follows the same API. Such a bespoke model is needed for the DRM.
• A MatrixModel needs to be written to handle DRM. Something like the following should be a good start:

  
  from astropy.modeling import Fittable1DModel, Parameter

class MatrixModel(Fittable1DModel): matrix = Parameter(description="The matrix with which to multiply the input.", fixed=True)

@staticmethod
def evaluate(x, matrix):
    return matrix @ x  # Requires input must have a specific dimensionality

* The input of the total model needs to be fixed to the spectral coords of the photon axis of the DRM.  This can be done with the ``astropy.modeling.fix_inputs` function.
* Physcial models in `sunkit-spex` needs to be converted to astropy models, or at least functions which can be wrapped in the `custom_model` decorator.
  * I am thinking in particular of the `thermal_emission` model which is currently a class.  But this should be fixable.

### `Spectrum` object
* Uncertainty objects need to be written in line with the `astropy.nddata.NDUncertainty` API.  These are needed for the following uncertainty types:
  * Poisson
  * etc.
* The fitter can check the type of this on the `Spectrum` object and use the correct error handling.
  * Perhaps some of the error handling methods can be attached directly onto these objects, similar to `StdDevUncertainty` and `VarianceUncertainty`.
  * A long term goal could be to upstream these to astropy, if we can make them follow the same API.

# Simultaneously fitting different models to multiple spectra OR same model to spectra on different spectral grids
* The `Fitter` should accept `tuples` of `(Spectrum, Model)` pairs. The fitting is then done by minimised the sum of the log-likelihoods of the different models to the different spectra.
Not sure whether this needs to be a separate collection fitter class, or if the same `Fitter` object should handle this case and the above one.

[DanRyanIrish]

The above is my suggested approach based on my current thinking. Open to changes as we discover more about implementing this in practice.

[DanRyanIrish]

After our productive discussion yesterday, I thought it would be helpful to clarify what I took from others' contributions, and how my thoughts have evolved on reflection. The main concerns I remember being raised were:

• Providing a tool that minimises the chance of user error;
• Linking the data and SRM as they are intrinsically related, and changing one (e.g. clipping spectral range, time range, etc.) will likely require a corresponding change to the other.
• Supporting a user-friendly way to make such changes, e.g. rebinning.

My feeling is that the reason we are finding these discussions so difficult, is that we are not sufficiently separating and abstracting the different functionalities we need and their requirements. Therefore, we often get trapped in discussions about a whole complex assortment of functionalities that the fitter should handle.

To make progress, I suggest that we focus only on the fundamentals of what we need for fitting, namely:

• A spectrum object
• A MatrixModel that follows the astropy modeling API and can therefore act as a component in a compound model (including physical models, optical path models, etc.)
  • This class acts as a wrapper around any object that supports the @ operator, i.e. matrix multiplication. This can be as simple as a numpy array, or a more complex instrument-specific SRM object with any number of methods, e.g. rebin.
• A fitter
  • This requires a spectrum object and a model that outputs a predicted spectrum on the same spectral grid, and the same dimensionality as the spectrum object. (fitter = Fitter((spectrum, model)))
  • No data manipulation is allowed. If new data or a new model is required, a new fitter instance must be generated.
  • However, masking certain ranges so they aren't included in the minimisation, or freezing/unfreezing, tying/untying parameters, etc. is supported.

This is the smallest set of infrastructure that enables users to do fitting.

If we want to enable a more user-friendly system that restricts chances for user error, I suggest we assign this to a future "fitting environment class", to be defined in the future. This could potentially take

• a fitter type
• a spectrum object
• a pre-srm model
• an srm
• a post-srm model
• sets of the above in cases of simultaneous fitting.

This object then support data/srm manipulation so long as the spectrum and srm support the same operations, e.g. .rebin(). Then, when it comes time to fit, a new fitter instance is generated with the self.spectrum object(s) and an on-the-fly-generated the model: model = self.pre_srm_model | self.srm | self.post_srm_model e.g. self.generate_fitter(*args) which doesn't something like self.fitter = self.fitter_type(self.spectrum, self.pre_srm_model | self.srm | self.post_srm_model) This could be done in a way in which various settings, e.g. which parameters are fixed, tied, etc., persist. Then the user could do fitter_env.fitter.fit() The API could also support multiple spectrum/model pairs for simultaneous fitting. The structure and API of this "fitting environment class" can be haggled over at a later stage. But crucially, this abstraction depends on the foundational infrastructure of

• spectrum
• model (including a matrix model)
• fitter

This should allow us to push ahead with the foundational infrastructure and deal with the "nice" user interface later.

[samaloney]

A MatrixModel that follows the astropy modeling API and can therefore act as a component in a compound model (including physical models, optical path models, etc.) This class acts as a wrapper around any object that supports the @ operator, i.e. matrix multiplication. This can be as simple as a numpy array, or a more complex instrument-specific SRM object with any number of methods, e.g. rebin.

Is very specific I would say need to abstract away to a Response that has some method forward, apply which is called or something e.g. no only / specifically matrix multiply

[DanRyanIrish]

I agree. There is probably a minimal API required by a representation of a response, e.g. the input and output spectral grid. But such a response object must be able to be a component in a compound astropy model, or be able to be wrapped in something that enables that.

My MatrixModel suggestion is simply a demonstration of how a response object could be wrapped into an astropy model component in the specific case where the response is represented by a matrix. But the fitting infrastructure should, in principle, support any astropy model.

[samaloney]

I guess the main point is that our Response object is more like a Spectrum object in terms of interface and function than a model. Also a large fraction of the time it's fixed/static so why are we trying to force into a model when are already have 90% of the required API and some code from the Spectrum/NDCube object?

We can always support wrapping a Response object with a model decorator if it needs to be varied in the fitting process but that not required right now to make progress.

[natsat0919]

A MatrixModel that follows the astropy modeling API and can therefore act as a component in a compound model (including physical models, optical path models, etc.) This class acts as a wrapper around any object that supports the @ operator, i.e. matrix multiplication. This can be as simple as a numpy array, or a more complex instrument-specific SRM object with any number of methods, e.g. rebin.

Is very specific I would say need to abstract away to a Response that has some method forward, apply which is called or something e.g. no only / specifically matrix multiply

I have seen the forward method in your toy-fitter code, what exactly is the purpose of the forward/apply methods?

[DanRyanIrish]

I guess the main point is that our Response object is more like a Spectrum object in terms of interface and function than a model. Also a large fraction of the time it's fixed/static so why are we trying to force into a model when are already have 90% of the required API and some code from the Spectrum/NDCube object?

I think this is an implementation detail. The epiphany I had after Wrocław is that philosophically, the response IS just a model component like any other, e.g. the thermal_emission model. It is a step in a chain of calculations that goes from photons emitted at the Sun, to counts measured by the detector. The fact that it may (or may not) be represented underneath as a matrix is a technicality, so long as it is possible to wrap the underlying matrix representation as a model, which the MatrixModel class proves it can be. This massively simplifies our design philosophy to 3 objects, i.e. spectrum + model + fitter.

Of course, we know that in practice, most users will want to use a response represented by a matrix. So we can provide convenience features for them to convert their matrix into a model, i.e. the MatrixModel class or an equivalent. We also know users may want to link the spectrum and response so that both are altered consistently, or forego the hassle of manually building the whole model themselves. Therefore, we can provide a sandbox in which these things can be done with reduced chance of user error and effort. But these are all convenience features built on top of a simpler, more powerful, albeit more error-prone, 3-object foundational API of spectrum + model + Fitter.