2D Rectangular Function Fitting Not Working

[William Steinberger]

Hello,

I am attempting to fit a data set with a 2D rectangular function. The data looks something like: 

import numpy as np
import matplotlib.pyplot as plt

from scipy.special import erf
import lmfit
from lmfit import Model

X,Y = np.meshgrid(x_Data,y_Data)
   

Max_Location_Y = np.where(z_Data==np.max(z_Data))[0][0]
Max_Location_X = np.where(z_Data==np.max(z_Data))[1][0]
   
X_3D = X[int(Max_Location_Y-10):int(Max_Location_Y+10)][:,int(Max_Location_X-10):int(Max_Location_X+10)]
Y_3D = Y[int(Max_Location_Y-10):int(Max_Location_Y+10)][:,int(Max_Location_X-10):int(Max_Location_X+10)]
Z_3D = z_Data[int(Max_Location_Y-10):int(Max_Location_Y+10)][:,int(Max_Location_X-10):int(Max_Location_X+10)]

fig = plt.figure(figsize = (10, 7))
ax = plt.axes(projection ="3d")
ax.scatter3D(X_3D, Y_3D, Z_3D/np.max(Z_3D), c='black', label="Data")
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Intensity')
plt.legend()
plt.show()
plt.close()

I get the following scatter plot:

I am then trying to fit the data to a 2D rectangular error function distribution: 

def Rectangle2D(x, y, amp, cen1x, cen2x, sig1x, sig2x ,cen1y, cen2y, sig1y, sig2y):
       
        Alpha_1x = (x-cen1x)/sig1x
        Alpha_2x = (x-cen2x)/sig2x
       
        Alpha_1y = (y-cen1y)/sig1y
        Alpha_2y = (y-cen2y)/sig2y
   
        Fit_Exp_Data = amp*(erf(Alpha_1x)-erf(Alpha_2x))*(erf(Alpha_1y)-erf(Alpha_2y))
       
        return Fit_Exp_Data.flatten()

gmodel = Model(Rectangle2D)
 result = gmodel.fit(Z_3D.flatten(), x=X_3D.flatten(), y=Y_3D.flatten(), amp=np.max(Z_3D), cen1x=168, cen2x=172, sig1x=1.0, sig2x=1.0, cen1y=218, cen2y=222, sig1y=1.0, sig2y=1.0)
lmfit.report_fit(result)

I then get the following error: 

"""

    result = gmodel.fit(Z_3D.flatten(), x=X_3D.flatten(), y=Y_3D.flatten(), amp=np.max(Z_3D), cen1x=168, cen2x=172, sig1x=1.0, sig2x=1.0, cen1y=218, cen2y=222, sig1y=1.0, sig2y=1.0)
  File "C:\Users\Z004C8MS\Anaconda3\lib\site-packages\lmfit\model.py", line 984, in fit
    params[name].set(value=p)
  File "C:\Users\Z004C8MS\Anaconda3\lib\site-packages\lmfit\parameter.py", line 668, in set
    self.value = value
  File "C:\Users\Z004C8MS\Anaconda3\lib\site-packages\lmfit\parameter.py", line 824, in value
    if self._val > self.max:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

"""

I have successfully done this with a 2D Gaussian:

model = lmfit.models.Gaussian2dModel()
params = model.guess(z_Data.flatten(), X.flatten(), Y.flatten())  # method cannot deal with 2D data

result = model.fit(z_Data, x=X, y=Y, params=params)
    lmfit.report_fit(result)

I cannot seem to implement the same thing for the lmfit Model for a 2D rectangle. Any help or suggestions would be greatly appreciated! 

Very Respectfully, 

Will 

[Matt Newville]

Hi William,

On Mon, Aug 29, 2022 at 1:21 PM William Steinberger <wm...@umich.edu> wrote:

Hello,

I am attempting to fit a data set with a 2D rectangular function. The data looks something like: 

One of the reasons that we continually ask for a minimal, complete example that demonstrates a problem is that we really do not have the time to try to fully understand your code.  You did not include complete, runnable code.  Specifically, `x_Data`, and `y_Data` are undefined so your code would not run past the first few lines.  

The error message from `parameter` is probably coming because you have defined a parameter that has a non-scaler value.  I would guess that is from `amp=np.max(Z_3D)`.  This comes from `z_Data` which is also undefined.   Maybe that's the problem.   Without runnable code, we can only guess.

I strongly recommend a) simplifying the code and making it more readable, and b) making a `Parameters` objects and evaluating the model function before trying a fit. 

--Matt