Spline approximation with lmfit

[Azat Khadiev]

Hello all,

I suppose that somebody has already faced the issue of spline approximation of the curve with lmfit package. I would like to use the spline function in my composite model to approximate the background in diffractogram.  I have found one class in scipy (Bspline), that could be useful for this purpose, but it requires knots and spline coefficients parameters as lists and it is not clear for me how to use lists as parameters for lmfit. I have found that M.Newville has some example code (https://groups.google.com/g/lmfit-py/c/tuSjIFUyLrg/m/NJPIHdeSAAAJ) for spline approximation, I would be appreciated if he could share it with me.

Regards, Azat. 

[Matt Newville]

On Thu, Apr 15, 2021 at 5:08 PM Azat Khadiev <azat.k...@gmail.com> wrote:

Hello all,

I suppose that somebody has already faced the issue of spline approximation of the curve with lmfit package. I would like to use the spline function in my composite model to approximate the background in diffractogram.  I have found one class in scipy (Bspline), that could be useful for this purpose, but it requires knots and spline coefficients parameters as lists and it is not clear for me how to use lists as parameters for lmfit. I have found that M.Newville has some example code (https://groups.google.com/g/lmfit-py/c/tuSjIFUyLrg/m/NJPIHdeSAAAJ) for spline approximation, I would be appreciated if he could share it with me.

Here's a simple example using just `scipy.interpolate.splev`.  Since you know at least some about splines, I'll skip much documentation. It's a pretty straightforward use of `splev` and B-splines (here, cubic).  To be clear, you could use `splrep` to generate a spline representation, which does handle endpoints well.  

Here, I just sort of hardwire in the offsets and the extra knots and so on.  You can explore that and play with the spline functions all you like. The point here is that you want the spline coefficients to be the variables. Also note that the values are such that you also can get decent guesses from the data to be interpolated (assuming you have that).

#!/usr/bin/env python

import numpy as np
from scipy.interpolate import splev
from lmfit import Parameters, minimize, fit_report
import matplotlib.pyplot as plt

FMT_COEF = 'coef_%2.2i'

def index_of(arrval, value):
    """return index of array *at or below* value """
    if value < min(arrval):
        return 0
    return max(np.where(arrval <= value)[0])

def spline_eval(params, xdata, xknots, degree=3):
    coefs = [params[FMT_COEF % i].value for i in range(ncoefs-degree-1)]
    return splev(xdata, [xknots, coefs, degree])

def residual(params, xdata, ydata, xknots, degree=3):
    return ydata - spline_eval(params, xdata, xknots, degree=degree)

xdata = np.linspace(0, 100, 201)
ydata = 8 + xdata/2.0 + 2*np.sin(xdata/2.5) + 5 * np.sin((xdata-2)/8)
ydata += np.random.normal(scale=0.1, size=len(xdata))

# create a limited number of spline "knots": 
#  they don't need to be evenly spaced in "x", but they are here for simplicity
ncoefs = 12
xknots = np.linspace(0, 100, ncoefs)

# now build ncoefs-4 (for 3rd-degree spline) parameters.
# note that you can get a pretty decent guess of the value
# by using the data point value close to x=xknot[i],
params = Parameters()
for i in range(ncoefs-4):
    params.add(name=FMT_COEF % i, value=ydata[index_of(xdata, xknots[i])])

print(params)
yinit = spline_eval(params, xdata, xknots)

result = minimize(residual, params, kws=dict(xdata=xdata, ydata=ydata, xknots=xknots))

print(fit_report(result, min_correl=0.5))
plt.plot(xdata, ydata, label='data')
plt.plot(xdata, yinit, label='init')
plt.plot(xdata, spline_eval(result.params, xdata, xknots), label='spline fit')
plt.legend()
plt.show()

#######

--Matt

[Azat Khadiev]

Thank you very much! I will try it!

суббота, 17 апреля 2021 г. в 21:13:34 UTC+2, Matt Newville:

[Azat Khadiev]

Dear Matt Newville,

I have tried this example code, it works well with the minimize() function. Thank you very much. In our case, I would like to use the spline function as a part of a composite function made of different Models (several Gaussians, etc). Since the number of coefficients and thus spline_eval() function's parameters can vary depending on the number of knots, it is not straightforward to use Model() class to wrap spline_eval() or splev() function. Does one need to explicitly define the number of function arguments, like:

def bspline(xdata, xknots, degree, c00, c01, c02, c03, c04):

   return splev(xdata, [xknots, [c00, c01, c02, c03, c04], degree])

every time we change the number of knots or degree to make it wrappable with Model() class?

Or there is some delicate way of doing that?

Regard, Azat.

понедельник, 19 апреля 2021 г. в 10:45:52 UTC+2, Azat Khadiev:

[Matt Newville]

Hi Azat, 

On Tue, Apr 20, 2021 at 6:01 AM Azat Khadiev <azat.k...@gmail.com> wrote:

Dear Matt Newville,

I have tried this example code, it works well with the minimize() function. Thank you very much. In our case, I would like to use the spline function as a part of a composite function made of different Models (several Gaussians, etc). Since the number of coefficients and thus spline_eval() function's parameters can vary depending on the number of knots, it is not straightforward to use Model() class to wrap spline_eval() or splev() function.

You're right.  In fact, as I was extracting that example from my "real-world" code, I also thought "hm, it would be nice if this could be a Model that handled all the ugly bits" and also "that could be really helpful for removing backgrounds in our integrated 1D XRD patterns".   So, maybe we're on the same page and this is something that we can add.

Does one need to explicitly define the number of function arguments, like:

def bspline(xdata, xknots, degree, c00, c01, c02, c03, c04):

   return splev(xdata, [xknots, [c00, c01, c02, c03, c04], degree])

every time we change the number of knots or degree to make it wrappable with Model() class?

Sort of.  One possible example to follow would be the PolynomialModel (https://github.com/lmfit/lmfit-py/blob/master/lmfit/models.py#L280).   That allows a Polynomial to be defined with up to 8 variables, and sort of hides much of the ugly parts from the user.  A spline could be similar (but maybe with a larger upper limit - even perhaps 50).   Still, the x-points for the knots would need to be specified and would have to agree with the degrees of freedom in the spline.

Or there is some delicate way of doing that?

There have been a few requests to be able to build a Model that does not use the explicit arguments of the wrapped model function but instead use some other packing of the actual parameters like a "list of Parameter values" so that the model function would be like a "var-args" function.  That would take some real effort, but it might be worth considering. 

--Matt