dynamically removing parameters from an astropy.modeling.Fittable1DModel model

[moustakas]

I'm working on an emission-line fitting code and having some speed / timing issues using the astropy.modeling.Fittable1DModel Class. I would be grateful for any suggestions or thoughts from the group.

The figure below shows the crux of my problem / question.

The input optical spectrum is a pure emission-line spectrum of a star-forming galaxy at z=0.368 (the stellar continuum has been subtracted). From left to right, the three panels are zoomed into the (1) [OII] 3726,29 doublet; (2) the H-beta + [OIII] 4959,5007 lines; and (3) the H-alpha + [NII] 6548,84 doublet.

The data (and +/-1-sigma uncertainties) are shown as a shaded gray spectrum in every panel; the top row shows the best-fitting model (in orange) using EMLineModel1 (defined below) and the bottom row (in blue) shows the best-fitting model using EMLineModel2.

If inspired, you can reproduce this figure by grabbing the emlinefit.zip file, which has the data and code, and running it in a dedicated conda environment (and directory). Roughly:

conda create -y -n emlinefit python=3.9 astropy=4.2 matplotlib pyyaml
conda activate emlinefit
pip install seaborn
python emlinefit.py

EMLineModel1 and EMLineModel2 both define my astropy.modeling.Fittable1DModel, roughly following the instructions here. The only difference between the two is that EMLineModel1 has 40 emission lines[*] while EMLineModel2 has 24 emission lines.

However, the 16 emission lines in EMLineModel1 but not in EMLineModel2 are all outside of my observed-frame wavelength range (roughly 3600-9800 A) and therefore not constrained by the data. Nevertheless, the optimizer (astropy.modeling.fitting.LevMarLSQFitter) using EMLineModel1 takes nearly a minute longer to reach the same solution as EMLineModel2 (see the figure legend, below)!

My sample spans a range of redshift and so I don't know a priori which emission lines will be in the observed-frame wavelength range. (I have tried to fix the Parameters of the lines outside the wavelength range so they're not optimized. This change does lead to a significant speed-up using EMLineModel1---18 seconds vs 9.4 seconds using EMLineModel2---although this is still a factor of two slower.)

So, my question is: can I dynamically remove Parameters from EMLineModel1 based on the redshift of a given object so I can minimize the optimization time? In other words, can I transform EMLineModel1 (which has all the lines I could ever want to fit) into a lighter-weight EMLineModel2 (which only has the Parameters constrainable by the data) on the fly?

Thanks in advance for any ideas or suggestions (or corrections to my code!).

[*] Each emission line has three Parameters associated with it, an amplitude, a line-width (in km/s), and a velocity shift (in km/s) with respect to the nominal redshift. However, not all these parameters are independent. For example, the line-widths and velocity shifts of all the Balmer and helium lines are tied together; similarly, all the forbidden transitions are tied together, etc.

[github-actions]

Welcome to Astropy 👋 and thank you for your first issue!

A project member will respond to you as soon as possible; in the meantime, please double-check the guidelines for submitting issues and make sure you've provided the requested details.

If you feel that this issue has not been responded to in a timely manner, please leave a comment mentioning our software support engineer @embray, or send a message directly to the development mailing list. If the issue is urgent or sensitive in nature (e.g., a security vulnerability) please send an e-mail directly to the private e-mail feedback@astropy.org.

[moustakas]

For example, I was hoping something super-simple like the snippet below would work, but alas:

python remove-Parameter.py
Default MyModel: param1 attribute exists? - True
MyModel with param1 attribute deleted: param1 attribute exists? - True

I also checked if I could simply overwrite the param_names attribute of the model Class, but that also doesn't work. Fittable1DModel subclasses the astropy.modeling.Model class and I tried to track down in the source code how I could update all the attributes which depend on the number of parameters at instantiation (so I can dynamically shorten the set of parameters, based on the redshift of each object), but I couldn't quite figure it out.

I've also done some basic profiling of the fitting code and can post the results here if that would be helpful.

#!/usr/bin/env python
from astropy.modeling import Fittable1DModel, Parameter
class MyModel(Fittable1DModel):
    param1 = Parameter(name='param1', default=1.0)
    param2 = Parameter(name='param2', default=2.0)
    def __init__(self, param1=param1.default, param2=param2.default, tossparam1=False, **kwargs):
        super(MyModel, self).__init__(param1=param1, param2=param2, **kwargs)
        if tossparam1:
            delattr(self, 'param1')
    @staticmethod
    def evaluate(x, p1, p2):
        return p1 + p2*x
if __name__ == '__main__':
    allparams = MyModel(tossparam1=False)
    print('Default MyModel: param1 attribute exists? - {}'.format(hasattr(allparams, 'param1')))

    noparam1 = MyModel(tossparam1=True)
    print('MyModel with param1 attribute deleted: param1 attribute exists? - {}'.format(hasattr(noparam1, 'param1')))

[pllim]

Your use case aside, not sure if I understand the concept of this in principle... If you remove stddev from Gaussian, is it still technically a Gaussian?

Perhaps there are other ways you can find your redshift without going down the modeling rabbit hole like this. Have you asked the spectroscopy community?

[moustakas]

Thanks for the feedback, @plim. I want to clarify that the galaxy redshift is known.

Below is a simpler illustration of the issue I've encountered. Imagine an intrinsic galaxy spectrum with only two emission lines (of varying and independent amplitude, width, and kinematic velocity shift with respect to the laboratory rest wavelength) and nothing else.

I create an astropy Model with free parameters for Line 1 and Line 2. However, because of the redshift of the galaxy, imagine that my observed spectrum only covers the wavelengths overlapping Line 2---I have no knowledge or ability to constrain the parameters of Line 1. (Of course, in a different galaxy, my observed-frame spectrum may cover Line 1 and/or Line 2, depending on the redshift.)

Now imagine that fitting the spectrum below takes 60 seconds. However, if I remove Line 1 from my Model, the fitting only takes, say, 5 seconds! That's the mystery I'm ultimately trying to understand.

Or, put another way: Is there any way for me to dynamically (on-the-fly) remove Line 1 from Model 1 for this object (or, equivalently, have the optimizer totally ignore the free parameters associated with Line 1 since they cannot be constrained by the data)?

I am fitting a sample with >10^6 objects, so compute time is very important, hence my desire to shave as much time from the optimization as possible.

I hope this explanation makes more sense.

[pllim]

Have you tried using compound models instead? https://docs.astropy.org/en/latest/modeling/compound-models.html

Removing/adding model instead of parameter might be easier.

[moustakas]

@plim, thanks but I'm pretty certain my fitting problem is perfectly suited for a multi-Parameter model because there are astrophysical constraints between parameters that are standard priors (e.g., the [OIII] 4959,5007 doublet ratio is 3:1 due to the atomic transition probabilities).

astropy.modeling.fitting.LevMarLSQFitter uses scipy.optimize.leastsq which uses MINPACK, so I know there's going to be an irreducible fitting time. However, my position is that fixed Parameters, i.e., parameters which are not being optimized (in my case, parameters outside the observed-frame wavelength range), should not increase the compute time significantly.

I did some basic profiling of my code to try to understand this issue using cProfile and snakeviz (see below). The time spent on scipy.optimize._minpack._lmdif is fine but it's suspicious that nearly all this time is spent on astropy.modeling.fitting._fitter_to_model_params. Digging a bit more, I found this comment:

"TODO: Most of this code should be entirely rewritten; it should not be as inefficient as it is."
https://github.com/astropy/astropy/blob/main/astropy/modeling/fitting.py#L1596-L1598

So it looks like this is the reason why the compute time scales with the number of Parameters even if the additional parameters are fixed / do not need to be optimized: every optimization loop calls LevMarLSQFitter.objective_function which itself calls _fitter_to_model_params.

[keflavich]

I think the feature you're looking for is the ability to hold one parameter fixed while performing the fit.

e.g. mymodel.fixed['param1'] = True.

Have you tested performance using that attribute?

[moustakas]

@keflavich, thanks for the suggestion.

Yes, I am setting fixed=True for parameters which I know are outside my wavelength range and therefore cannot be optimized. However, I found many places in the code which loops on all the parameters of the model (irrespective of whether they're fixed or not), which is partially why I was finding the fitting time to scale with the number of parameters in the model.

Profiling ended up being a very fruitful exercise. I found significant bottlenecks in modeling.fitting._fitter_to_model_params as well as in all the bound properties of the modeling.core.Model class, such as tied, fixed, etc. These properties each call modeling.utils._ConstraintsDict, which loops on all the parameters in the param_names array every time they're called, and, in particular, many times inside the optimization loop.

By subclassing the modeling.fitter.LevMarLSQFitter fitter and trying to remove as many of these bottlenecks as possible, I've reduced the compute time (for my specific problem, with all 40 emission lines of interest) from 18.0 seconds with the 4.2.1 version of astropy to 1.55 seconds, a factor of 11.6 speed-up.

I'm an old convert from IDL, so my code is by no means the best way to do this, but I'm including it below in case others find it useful . Additional comments and suggestions are more than welcome!

And here's the full-spectrum fit, for completeness:

from astropy.modeling.fitting import LevMarLSQFitter
class MyLevMarLSQFitter(LevMarLSQFitter):
    def __init__(self, model):
        import operator        
        from functools import reduce

        self.has_tied = any(model.tied.values())
        self.has_fixed = any(model.fixed.values())
        self.has_bound = any(b != (None, None) for b in model.bounds.values())

        _, fit_param_indices = self.model_to_fit_params(model)

        self._fit_param_names = [model.param_names[iparam] for iparam in np.unique(fit_param_indices)]

        self._fit_istart = []
        self._fit_iend = []
        self._fit_slice = []
        self._fit_ibounds = []
        self._tied_slice = []
        self._tied_func = []

        offset = 0
        for name in self._fit_param_names:
            slice_ = model._param_metrics[name]['slice']
            shape = model._param_metrics[name]['shape']
            # This is determining which range of fps (the fitted parameters) maps
            # to parameters of the model
            size = reduce(operator.mul, shape, 1)
            self._fit_slice.append(slice_)
            self._fit_istart.append(offset)
            self._fit_iend.append(offset + size)

            imin, imax = model.bounds[name]
            self._fit_ibounds.append(((imin, imax) != (None, None), imin, imax))
            offset += size            

        if self.has_tied:
            self._tied_param_names = [_param_name for _param_name in model.param_names if model.tied[_param_name]]
            for name in self._tied_param_names:
                self._tied_slice.append(model._param_metrics[name]['slice'])
                self._tied_func.append(model.tied[name])
        
        super().__init__()

    def objective_function(self, fps, *args):

        model = args[0]
        weights = args[1]
        self.fitter_to_model_params(model, fps)

        meas = args[-1]
        if weights is None:
            return np.ravel(model(*args[2: -1]) - meas)
        else:
            return np.ravel(weights * (model(*args[2: -1]) - meas))

    def fitter_to_model_params(self, model, fps):
        """Constructs the full list of model parameters from the fitted and constrained
        parameters.

        """
        parameters = model.parameters

        if not (self.has_tied or self.has_fixed or self.has_bound):
            # We can just assign directly
            model.parameters = fps
            return

        for name, istart, iend, islice, (ibounds, imin, imax) in zip(
                        self._fit_param_names, self._fit_istart, self._fit_iend,
                        self._fit_slice, self._fit_ibounds):
            values = fps[istart:iend]

            # Check bounds constraints
            if ibounds:
                if imin is not None:
                    values = np.fmax(values, imin)
                if imax is not None:
                    values = np.fmin(values, imax)
                    
            parameters[islice] = values
            setattr(model, name, values)

        # Update model parameters before calling ``tied`` constraints.

        # This has to be done in a separate loop due to how tied parameters are
        # currently evaluated (the fitted parameters need to actually be *set* on
        # the model first, for use in evaluating the "tied" expression--it might be
        # better to change this at some point
        if self.has_tied:
            for name, islice, func in zip(
                    self._tied_param_names, self._tied_slice,
                    self._tied_func):
                #value = model.tied[name](model)

                # To handle multiple tied constraints, model parameters
                # need to be updated after each iteration.
                value = func(model)
                parameters[islice] = value
                setattr(model, name, value)
                
    def model_to_fit_params(self, model):
        """Convert a model instance's parameter array to an array that can be used with
        a fitter that doesn't natively support fixed or tied parameters.  In
        particular, it removes fixed/tied parameters from the parameter array.
        These may be a subset of the model parameters, if some of them are held
        constant or tied.

        """
        fitparam_indices = list(range(len(model.param_names)))
        if self.has_fixed or self.has_tied:
            params = list(model.parameters)
            param_metrics = model._param_metrics
            for idx, name in list(enumerate(model.param_names))[::-1]:
                if model.fixed[name] or model.tied[name]:
                    slice_ = param_metrics[name]['slice']
                    del params[slice_]
                    del fitparam_indices[idx]
            return (np.array(params), fitparam_indices)
        return (model.parameters, fitparam_indices)