Renamed parameters in Gaussian and added a cov_matrix kwarg - #1022

@eteq @nden

I like the name changes - more PEP8 compliant and easier to read.

Now in the Gaussian2D initializer: I have a suggestion for how to formulate this more cleanly - I see @adrn's concerns that there are a lot of ways to do the same thing. But you can make all the different ways of specifying the shape parameters be properties on the object. You can leave the internal representation alone (I think it's x_sigma, y_sigma, and theta, right?), and generate everything else on demand as properties (or set as properties). Then it's ok to have tons of different representations. __init__ can delegate to the property setters, and then instead of nested if/else statements, you can just enumerate all the valid combinations and fail once it reaches else.

As for the generic nd-gaussian, I suggest you implement that as a separate class (NDGaussian?) and a separate PR. There are plenty of people who just want a 1D gaussian, and lots of cases where you want the security of knowing it's 2D, as well (images). Once you get past that, n becomes arbitrary and it makes sense to have it be a general vector.

Also, there's a sphinx build failure - looks like you didn't update things that are in the docs for the new name.

I agree with @eteq that NDGaussian should be a separate model, there are plenty of use cases for simple 1D and 2D gaussians.

I have no objections to changing the parameter names in the gaussian models. I needed a test case for non-linear fitters, picked gaussian and looked up the names IRAF is using. Although 'center' is clearer to me than 'mu' I'll agree to calling it 'mu'.

You are right, no particular reason to keep 'ratio', an "irafism" again. But I am surprised 'fwhm' is not favored here as it seems to be the natural choice in some fields of astronomy. Anyway, I'll leave it to real users to decide.

No objection to the param_* changes either. This makes it a bit harder to work with in an interactive session but they are more descriptive this way.

The examples in the docs need to be changed as well.

Another reference point: gauss1d and gauss2d parameter names in Sherpa.

As a user I'm +1 on using sigma instead of fwhm (but provide the constant in astropy.constants to convert back and forth and maybe even mention that this constant is available in the docs.

And -1 on mu, I'd prefer pos as Sherpa does or center as @nden did. mu is just a greek letter that means different things in different contexts, whereas pos or center describe the parameter.

I don't think we need to choose which is the "best" between sigma or fwhm - both are common in various fields. That's why I said we should support all of them using properties and accept any valid combination in the initializer.

I'm +0 on center over mu - there I think most people will understand either way. (But I'm -1 on pos, which is more ambiguous.)

Personally, I prefer mean, but it doesn't make sense to me to have mean and sigma -- instead it should be mean and std_dev or something like that. So that's why I chose mu -- not sure why sigma is less opaque than mu, though, since it too is just a greek letter that means different things in different contexts :)

In fact, what I really prefer is mean and variance -- because if NDGaussian expects a covariance matrix, it is a more natural 1D extension of that.

@eteq I am not sure I understand your suggestion about handling the parameters through properties. If we take Gaussian1D as an example (ignoring parameter names now) , the relevant parameters are

 def __init__(self, amplitude, xcen, fwhm=None, xsigma=None)

I think @adrn would like to see them as:

 def __init__(self, amplitude, xcen, xsigma)

If I understand you correctly, you are proposing to keep the latter signature but allow users to get/set the fwhm using properties.
If so, I don't see how this makes things easier for users. Someone who estimates the fwhm visually would have to calculate the sigma to set the initial values. This may be a bit confusing. And someone who's used to supplying fwhm as a parameter (any IRAF user) may even get this wrong initially. The current parameters force users to be specific and so are less error prone.

Also, it may not be obvious, so I'll mention it here. All model parameters (those that are in the parnames list) are automatically turned into properties in the base Model class. This is done through core._parameterProperty.

I've been thinking more about this.
First, the behaviour of creating a Parameter from an array that @adrn reports is a bug, I'll fix it.
But cov_matrix does not need to be of type Parameter.

My intention from the start was to be able to create models for which one or more parameters may be of type array. However it is not clear to me what it means to fit a parameter which is an array. What I had in mind is something similar to what @adrn has done with cov_matrix. Such parameters are exposed at the user level for convenience but the model class extracts from them meaningful parameters of type number to be used for fitting.

With that in mind and this Gauss2D definition as an example, here's how parameters work.

class Gauss2D():
    param_names = ['amplitude', 'x_mu', 'y_mu', 'x_sigma', 'y_sigma', 'theta']

    def __init__(self, 'amplitude', x_mu, y_mu, x_fwhm=None, y_fwhm=None,
                     x_sigma=None, y_sigma=None, theta=0.0, cov_matrix=None)

param_names defines the parameters which the model uses internally (for fitting and evaluation). These must be turned to type Parameter in __init__. Also _ParameterProperty creates automatically properties for them. (Side note: the setter is especially important because it validates each parameter.)
With the above definition it is possible to do:

Gauss2D(amplitude=10, x_mu=3, y_mu=3, x_sigma=.4, y_sigma=.3, theta=13)

or

Gauss2D(amplitude=10, x_mu=3, y_mu=3, x_fwhm=.4, y_fwhm=.3, theta=13)

or

Gauss2D(amplitude=10, x_mu=3, y_mu=3, cov_matrix=my_matrix, theta=13)

Gauss2D will have properties (created manually, not through _ParameterProperty) for cov_matrix, x_fwhm and y_fwhm and the setters MUST create x_sigma and y_sigma as type Parameter. I guess this is what @eteq was describing except I thought @adrn did not like fwhm as an optionl input.

How does this sound? Is it flexible enough?

@nden - yes, that's what I was thinking (although I think the third one shouldn't have theta, right? That is, theta is computed from the off-diagonals of cov_matrix?)

I think I was being confusing because I was mostly talking about the properties, rather than the initializer. I was saying that we should accept all of the options as possible initializer arguments (as you are showing above), but they should also all be accessible (or settable) as properties. But I agree that that any non-overdetermined combination should be acceptable in the initializer (including e.g. x_sigma and y_fwhm as long as y_sigma and x_fwhm are not given).

By the way, @adrn, once the initializer plan is finalized, given that you'll be changing the docstrings anyway, can you also add to the documentation which direction theta goes? There is a long and sordid history of erroneous science due to confusing the direction of position angles!

These changes touch every file in the package, we should merge this PR first to minimize future work.
Let's merge this soon and move on.

So to sum up, I think there's agreement on changing these names:

• Gauss1DModel --> Gaussian1DModel
• Gauss2DModel --> Gaussian2DModel
• psets --> param_sets
• paramdim --> param_dim
• parnames - param_names

The Gaussian2DModel initalizer is now:

class Gaussian2DModel():
    param_names = ['amplitude', 'x_mu', 'y_mu', 'x_sigma', 'y_sigma', 'theta']

    def __init__(self, 'amplitude', x_mu, y_mu, x_fwhm=None, y_fwhm=None,
                     x_sigma=None, y_sigma=None, theta=0.0, cov_matrix=None)

Internally the model works with the parameters listed in param_names and properties are created for them automatically. There needs to be a property for each of x_fwhm, y_fwhm, and cov_matrix. These properties set the internal parameters.

Is this acceptable to everyone?

And finally, two options for the parameter names emerged:

[amplitude, center, sigma]

or

[amplitude, mean, variance]

Either of these looks good to me. Any final votes?
If there are no strong feelings either way, I'd say @adrn should pick one and go with it.

What needs to be done before this is merged is minor fixes to the documentation and adding the properties.

I think using variance = width ** 2 instead of sigma = width as a parameter is bad.

It's a non-linear relationship variance = sigma ** 2 and while best-fit values are simply related by that formula, errors on best-fit values are not and can be very tricky near the hard limit sigma = 0.

My argument for using sigma and not variance is that it is more common to want a best-fit value and error for the width instead of the width ** 2 so that's e.g. what Sherpa and the other modelling packages I know do.

I have a strong preference for sigma over variance as that is what 99% of gaussian model fits I've seen before use.

Ok, fine by me. But I don't think mean and sigma make sense together -- I think either [mu, sigma] or [mean, stddev]

@adr I think I argued against mu before ... but actually either mu, sigma or mean, stddev would be fine with me ... go ahead with what you like better.

Travis failure seems unrelated...

Just restarted the travis failed build - probably that will fix it (as @adrn said, it looked like just a timeout, nothing real).

Assuming it passes after that, this looks good to merge to me.

Cool, build passed!

Thanks @adrn and @eteq ! I'll merge this.