numpysympylmfit

numpy ndarray error in lmfit when mdel is passed using sympy


I got the following error:

<lambdifygenerated-1>:2: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.return numpy.array((A1exp(-1/2(x - xc1)**2/sigma1**2), 0, 0))

Here I have just one model but this code is written for model combination in fitting by the lmfit Please kindly let me know about it.

import matplotlib.pyplot as plt
import numpy as np
import sympy
from sympy.parsing import sympy_parser
import lmfit

gauss_peak1 = sympy_parser.parse_expr('A1*exp(-(x-xc1)**2/(2*sigma1**2))')
gauss_peak2 = 0
exp_back = 0

model_list = sympy.Array((gauss_peak1, gauss_peak2, exp_back))
model = sum(model_list)
print(model)



model_list_func = sympy.lambdify(list(model_list.free_symbols), model_list)
model_func = sympy.lambdify(list(model.free_symbols), model)


np.random.seed(1)
x = np.linspace(0, 10, 40)
param_values = dict(x=x, A1=2, sigma1=1, xc1=2)
y = model_func(**param_values)
yi = model_list_func(**param_values)
yn = y + np.random.randn(y.size)*0.4

plt.plot(x, yn, 'o')
plt.plot(x, y)


lm_mod = lmfit.Model(model_func, independent_vars=('x'))
res = lm_mod.fit(data=yn, **param_values)
res.plot_fit()
plt.plot(x, y, label='true')
plt.legend()
plt.show()

Solution

  • lmfit.Model takes a model function that is a Python function. It parses the function arguments and expects those to be the Parameters for the model.

    I don't think using sympy-created functions will do that. Do you need to use sympy here? I don't see why. The usage here seems designed to make the code more complex than it needs to be. It seems you want to make a model with a Gaussian-like peak, and a constant(?) background. If so, why not do

    from lmfit.Models import GaussianModel, ConstantModel
    
    model = GaussianModel(prefix='p1_') + ConstantModel()
    params = model.make_params(p1_amplitude=2, p1_center=2, p1_sigma=1, c=0)
    

    That just seems way easier to me, and it is very easy to add a second Gaussian peak to that model.

    But even if you have your own preferred mathematical expression, don't use that as a sympy string, use it as Python code:

    def myfunction(x, A1, xc1, sigma1):
        return A1*exp(-(x-xc1)**2/(2*sigma1**2))
    

    and then

    from lmfit import Model
    mymodel = Model(myfunction)
    params = mymodel.guess(A1=2, xc1=2, sigma1=1)
    

    In short: sympy is an amazing tool, but lmfit does not use it.