gaussian2d doesn't actually give a 2d-matrix

[bguttel]

It seems to me that the gaussian2d function under lineshapes just multiply two 1D-gaussians, and since this is the normal numpy multiplication it is only element-wise and hence the result is also one-dimensional (with the same shape).
In practice it is the diagonal of the expected matrix.

I would say the multiplication needs to be the numpy.outer multiplication. If that's true I can try and write it down..

[newville]

@bguttel is https://lmfit.github.io/lmfit-py/examples/example_two_dimensional_peak.html helpful?

[bguttel]

@bguttel is https://lmfit.github.io/lmfit-py/examples/example_two_dimensional_peak.html helpful?

Actually this example is very confusing, because its author used the scipy.interpolate.griddata function. So most of the plot is just linear interpolation of the original data and its fit.
The reason it looks so convincing is because the "2d-gaussian" is calculated on plenty of random points of the grid.

Edit: in this example x,y and z all have the same shape (10000,) which of course doesn't make sense.

[newville]

@bguttel Right, I see your point that requiring x and y to be the same size is not good. This function probably ought to return a 2D array. @reneeotten might have some thoughts on this too, but maybe we can settle on a more "obvious" 2D function form and include another example or two, say lorentzian2d?

[bguttel]

Definitely.
I encountered this problem when I was trying to fit actual 2d-data using this built-in model, which failed. I think the function should look similar to this:
def matrix_gaussian2d(x, y, amplitude=1.0, centerx=0.0, centery=0.0, sigmax=1.0, sigmay=1.0): """Return a flatten 2-dimensional Gaussian function.... """ z = amplitude*(np.outer(gaussian(y, amplitude=1, center=centery, sigma=sigmay), gaussian(x, amplitude=1, center=centerx, sigma=sigmax))) return z.flatten()

It will be nice to return (and receive as data-to-be-fitted) a non-flattened 2d-data, but I am not sure lmfit supports that?..

[newville]

@bguttel I think it is reasonable to return a 2D as you suggest and then make sure that any 2D (or other) model gets flattened in the fit -- that should happen already.

[reneeotten]

I thought it was returning a 2D array, but a cursory look indeed suggests that's not happening right now... I can certainly take a look, but perhaps @mpmdean who added this can elaborate the current implementation?

[mpmdean]

gaussian2d does give you a 2D array, but only if you pass it a 2D array.

import numpy as np
from lmfit.lineshapes import gaussian2d

X, Y = np.meshgrid(np.linspace(-1, 1, 30),
                   np.linspace(-1, 1, 30))

Z = gaussian2d(X, Y)

Z.shape is then (30, 30)

The advantage of this is that you can fit data even if it is not on a regular grid. If the data is not on a regular grid one cannot represent it in terms of the two 1D arrays of length N and M to make up a 2D NxM array.

I'm open to input, but my understanding is that scripts with the meshgrid step will all work perfectly?

[newville]

@mpmdean thanks for the clarification - that does seem like a compelling use case. I think perhaps the only example we have is slightly confusing. Perhaps the need to use numpy.meshgrid to turn "simple 1D x and y arrays" into the proper inputs should be clarified in the example?

[mpmdean]

I would agree that the example would be improved by addressing this more explicitly. Now is not a good time for me in terms of availability to address this, I'm afriad.

[newville]

@mpmdean oh no worries -- and I'm not able to do it any time soon either! Maybe @bguttel would be willing to provide or improve the example?

[bguttel]

Sorry for the late response. Following @mpmdean explanation, I think works it would be helpful to add a semi-example with 2d data (Z) and x and y 1d vectors as follows:

X, Y = np.meshgrid(x, y)
model = lmfit.models.Gaussian2dModel()
params = model.guess(Z.flatten(), X.flatten(), Y.flatten())
result = model.fit(Z, x=X, y=Y, params=params)

Notice that I also had to flatten the matrices since the guess_from_peak2d function only works that way.