0

I am trying to fit some experimental data (x and y) with a custom function (Srt) and using scipy.optimize.curve_fit():

Reading the data and defining the function, using dummy values (10,10) for Km and Vmax (which are to be determined using the curve fit) works fine, as long as I use np.asarray():

from scipy.special import lambertw
from scipy.optimize import curve_fit
import numpy as np
import scipy

def Srt(t,s,Km,Vmax):
    print("t",type(t))
    print("t",t)
    print("last element of t:",t[-1])
    print("s",type(s))
    print("s",s)
    print("last element of s:",s[-1])
    Smax = s[-1] # Substrate concentration at end of reaction 
    t0 = t[0]    # time=0 (beginning of reaction)
    s0 = s[0]    # Substrate concentration at time = 0 (beginning of reaction)
    E = np.exp(((Smax - s0) - Vmax*(t+t0))/Km)
    L = lambertw(((Smax - s0)/Km)*E)
    y = Smax - Km*L
    return y

x=[2.780000e-03,2.778000e-02,5.278000e-02,7.778000e-02,1.027800e-01
,1.277800e-01,1.527800e-01,1.777800e-01,2.027800e-01,2.277800e-01
,2.527800e-01,2.777800e-01,3.027800e-01,3.277800e-01,3.527800e-01]

y=[0.44236,0.4308,0.42299,0.41427,0.40548,0.39908,0.39039,0.3845,0.37882
,0.37411,0.36759,0.36434,0.35864,0.35508,0.35138]

xdata = np.asarray(x)
ydata = np.asarray(y)
Srt(xdata, ydata,10,10)

If I do not use np.asarray, I get a "Type Error":

Srt(x, y,10,10)

TypeError

When I continue to use curve_fit to make the fit for Vmax and Km with:

parameters, covariance = scipy.optimize.curve_fit(Srt, xdata, ydata)

I get into trouble:

If I understand the error message correctly, for some reason, the array ydata is not an array anymore when it is read in as s ?!?

ydata_S_Error

What do I have to change in my code so that I can work with my function Srt and curve_fit ?

EDIT: Full output of code:

t <class 'numpy.ndarray'>
t [0.00278 0.02778 0.05278 0.07778 0.10278 0.12778 0.15278 0.17778 0.20278
 0.22778 0.25278 0.27778 0.30278 0.32778 0.35278]
last element of t: 0.35278
s <class 'numpy.ndarray'>
s [0.44236 0.4308  0.42299 0.41427 0.40548 0.39908 0.39039 0.3845  0.37882
 0.37411 0.36759 0.36434 0.35864 0.35508 0.35138]
last element of s: 0.35138
t <class 'numpy.ndarray'>
t [0.00278 0.02778 0.05278 0.07778 0.10278 0.12778 0.15278 0.17778 0.20278
 0.22778 0.25278 0.27778 0.30278 0.32778 0.35278]
last element of t: 0.35278
s <class 'numpy.float64'>
s 1.0

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-23-5ce34d06b849> in <module>
     33 #then the problems start
     34 
---> 35 parameters, covariance = scipy.optimize.curve_fit(Srt, xdata, ydata)

~\anaconda3\lib\site-packages\scipy\optimize\minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    761         # Remove full_output from kwargs, otherwise we're passing it in twice.
    762         return_full = kwargs.pop('full_output', False)
--> 763         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    764         popt, pcov, infodict, errmsg, ier = res
    765         ysize = len(infodict['fvec'])

~\anaconda3\lib\site-packages\scipy\optimize\minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    386     if not isinstance(args, tuple):
    387         args = (args,)
--> 388     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    389     m = shape[0]
    390 

~\anaconda3\lib\site-packages\scipy\optimize\minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     24 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     25                 output_shape=None):
---> 26     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     27     if (output_shape is not None) and (shape(res) != output_shape):
     28         if (output_shape[0] != 1):

~\anaconda3\lib\site-packages\scipy\optimize\minpack.py in func_wrapped(params)
    461     if transform is None:
    462         def func_wrapped(params):
--> 463             return func(xdata, *params) - ydata
    464     elif transform.ndim == 1:
    465         def func_wrapped(params):

<ipython-input-23-5ce34d06b849> in Srt(t, s, Km, Vmax)
     10     print("s",type(s))
     11     print("s",s)
---> 12     print("last element of s:",s[-1])
     13     Smax = s[-1] # Substrate concentration at end of reaction
     14     t0 = t[0]    # time=0 (beginning of reaction)

IndexError: invalid index to scalar variable.

EDIT 2 FULLY Functional code, thanks to Jonathan Weine. Fit is suboptimal due to "bad" experimental data, I am playing around with my full dataset now :D

from scipy.special import lambertw
from scipy.optimize import curve_fit
import numpy as np
import scipy

def Srt(t, Smax: float, s0: float, Km: float, Vmax: float):
    t0 = t[0]    # time=0 (beginning of reaction)
    E = np.exp(((Smax - s0) - Vmax*(t+t0))/Km)
    # L = lambertw(((Smax - s0)/Km)*E)  # this apparently can be complex which causes another Error
    L = np.abs(lambertw(((Smax - s0)/Km)*E))
    y = Smax + Km*L
    return y

y = [0.44236,0.4308,0.42299,0.41427,0.40548,0.39908,0.39039,0.3845,0.37882
,0.37411,0.36759,0.36434,0.35864,0.35508,0.35138,0.34748,0.34437,0.34143
,0.3391,0.3348,0.33345,0.31404,0.30212,0.29043,0.28026,0.27331,0.26672
,0.26187,0.25645,0.25208,0.24736,0.244,0.24056,0.23798,0.23359,0.23138
,0.22845,0.22566,0.22384,0.22112,0.21894,0.21672,0.21466,0.21316,0.21209
,0.20941,0.20823,0.20687,0.2056,0.20324,0.20266,0.20095,0.19935,0.19895
,0.19746,0.19616,0.19486,0.19419,0.19382,0.19301,0.19085,0.19108,0.19024
,0.18933,0.18839,0.18706,0.18643,0.18623,0.18569,0.18469,0.18381,0.18341
,0.18331,0.18324,0.18222,0.18106,0.18039,0.18022,0.17906,0.17935,0.17842
,0.17834,0.1781,0.17731,0.17704,0.1766,0.17654,0.1761,0.17568,0.1744
,0.17453,0.17393,0.17325,0.17329,0.17302,0.17347,0.17344,0.17233,0.17228
,0.17208,0.17177,0.1712,0.17076,0.171,0.17043,0.17057,0.17003,0.16965
,0.16923,0.16944,0.16898,0.16879,0.16809,0.16821,0.16794,0.16831,0.16779
,0.16805,0.16765,0.16762,0.16695,0.16694,0.1669,0.16642,0.16583,0.166
,0.16625,0.16575,0.1658,0.16553,0.16565,0.1654,0.16419,0.16487,0.16467
,0.16452,0.16433,0.16468,0.16423,0.16427,0.16372,0.16388,0.16388,0.16394
,0.16382,0.1631,0.16353,0.1638,0.16304,0.163,0.16296,0.16295,0.16284
,0.16275,0.16214,0.16243,0.16211,0.16207,0.16185,0.16187,0.16176,0.16168
,0.16195,0.16138,0.16177,0.16121,0.16163,0.16121,0.161,0.16114,0.16122
,0.16096,0.16105,0.16102,0.16068,0.16031,0.16028,0.16051,0.16045,0.16017
,0.15977,0.15927,0.16007,0.15953,0.15933,0.1596,0.15911,0.15903,0.15884
,0.15856,0.15889,0.15888,0.15861,0.15849,0.158,0.15822,0.15776,0.15759
,0.15734,0.15757,0.15718,0.15699,0.15747,0.15692,0.15701,0.15715,0.15675
,0.15732,0.15687,0.15659,0.15664,0.15635,0.15633,0.15591] 
#csvFile.iloc[0:500,9] 

x = [2.780000e-03,2.778000e-02,5.278000e-02,7.778000e-02,1.027800e-01
,1.277800e-01,1.527800e-01,1.777800e-01,2.027800e-01,2.277800e-01
,2.527800e-01,2.777800e-01,3.027800e-01,3.277800e-01,3.527800e-01
,3.777800e-01,4.027800e-01,4.277800e-01,4.527800e-01,4.777800e-01
,5.027800e-01,7.538900e-01,1.003890e+00,1.253890e+00,1.503890e+00
,1.753890e+00,2.003890e+00,2.253890e+00,2.503890e+00,2.753890e+00
,3.003890e+00,3.253890e+00,3.503890e+00,3.753890e+00,4.003890e+00
,4.253890e+00,4.503890e+00,4.753890e+00,5.003890e+00,5.253890e+00
,5.503890e+00,5.753890e+00,6.003890e+00,6.253890e+00,6.503890e+00
,6.753890e+00,7.003890e+00,7.253890e+00,7.503890e+00,7.753890e+00
,8.003890e+00,8.253890e+00,8.503890e+00,8.753890e+00,9.003890e+00
,9.253890e+00,9.503890e+00,9.753890e+00,1.000389e+01,1.025389e+01
,1.050389e+01,1.075389e+01,1.100389e+01,1.125389e+01,1.150389e+01
,1.175389e+01,1.200389e+01,1.225389e+01,1.250389e+01,1.275389e+01
,1.300389e+01,1.325389e+01,1.350389e+01,1.375389e+01,1.400389e+01
,1.425389e+01,1.450389e+01,1.475389e+01,1.500389e+01,1.525389e+01
,1.550389e+01,1.575389e+01,1.600389e+01,1.625389e+01,1.650389e+01
,1.675389e+01,1.700389e+01,1.725389e+01,1.750389e+01,1.775389e+01
,1.800389e+01,1.825389e+01,1.850389e+01,1.875389e+01,1.900389e+01
,1.925389e+01,1.950389e+01,1.975389e+01,2.000389e+01,2.025389e+01
,2.050389e+01,2.075389e+01,2.100389e+01,2.125389e+01,2.150389e+01
,2.175389e+01,2.200389e+01,2.225389e+01,2.250389e+01,2.275389e+01
,2.300389e+01,2.325389e+01,2.350389e+01,2.375389e+01,2.400389e+01
,2.425389e+01,2.450389e+01,2.475389e+01,2.500389e+01,2.525389e+01
,2.550389e+01,2.575389e+01,2.600389e+01,2.625389e+01,2.650389e+01
,2.675389e+01,2.700389e+01,2.725389e+01,2.750389e+01,2.775389e+01
,2.800389e+01,2.825389e+01,2.850389e+01,2.875389e+01,2.900389e+01
,2.925389e+01,2.950389e+01,2.975389e+01,3.000389e+01,3.025389e+01
,3.050389e+01,3.075389e+01,3.100389e+01,3.125389e+01,3.150389e+01
,3.175389e+01,3.200389e+01,3.225389e+01,3.250389e+01,3.275389e+01
,3.300389e+01,3.325389e+01,3.350389e+01,3.375389e+01,3.400389e+01
,3.425389e+01,3.450389e+01,3.475389e+01,3.500389e+01,3.525389e+01
,3.550389e+01,3.575389e+01,3.600389e+01,3.625389e+01,3.650389e+01
,3.675389e+01,3.700389e+01,3.725389e+01,3.750389e+01,3.775389e+01
,3.800389e+01,3.825389e+01,3.850389e+01,3.875389e+01,3.900389e+01
,3.925389e+01,3.950389e+01,3.975389e+01,4.000389e+01,4.025389e+01
,4.050389e+01,4.075389e+01,4.100389e+01,4.125389e+01,4.150389e+01
,4.175389e+01,4.200389e+01,4.225389e+01,4.250389e+01,4.275389e+01
,4.300389e+01,4.325389e+01,4.350389e+01,4.375389e+01,4.400389e+01
,4.425389e+01,4.450389e+01,4.475389e+01,4.500389e+01,4.525389e+01
,4.550389e+01,4.575389e+01,4.600389e+01,4.625389e+01,4.650389e+01
,4.675389e+01,4.700389e+01,4.725389e+01,4.750389e+01,4.775389e+01
,4.800389e+01,4.825389e+01,4.850389e+01,4.875389e+01] 
# csvFile.iloc[0:500,0] 

xdata = np.array(x)
ydata = np.array(y)

parameters, covariance = scipy.optimize.curve_fit(Srt, xdata, ydata)
Km = parameters[1]
Vmax = parameters[0]
fit_y = Srt(xdata, ydata[-1],ydata[0], Km, Vmax)

print("Km: ", parameters[1], "Vmax: ", parameters[0])
plt.plot(xdata, fit_y, '-', color="green",linewidth=1)
plt.plot(xdata, ydata, 'o', color="red")
Improve this question
2
  • 1
    With numpy, you are using broadcasting numpy.org/doc/stable/user/basics.broadcasting.html, Without numpy, your calculation is float + list, which throws the error about the failed concatenation. Commented Dec 14, 2021 at 15:27
  • Post full tracebacks as text, not images. Your function should be callable with a scalar. Remove the inappropriate printouts, or change them to accept a scalar properly. Commented Dec 14, 2021 at 15:27

2 Answers 2

Reset to default
1

Please have a closer look at the documentation of the curve_fit function. Where it states that ydata must nominaly be the result of func(xdata... ). So the ydata that you hand to curve_fit is never passed as argument of the call of Srt as you indicated in the manual call.

Furthermore, the parameters to be estimated must have the same shape, which means that you have to define Smax and s0 as float input. I modified your example such that it actually runs:

from scipy.special import lambertw
from scipy.optimize import curve_fit
import numpy as np
import scipy

def Srt(t, Smax: float, s0: float, Km: float, Vmax: float):
    t0 = t[0]    # time=0 (beginning of reaction)
    E = np.exp(((Smax - s0) - Vmax*(t+t0))/Km)
    # L = lambertw(((Smax - s0)/Km)*E)  # this apparently can be complex which causes another Error
    L = np.abs(lambertw(((Smax - s0)/Km)*E))
    y = Smax - Km*L
    return y

x=[2.780000e-03,2.778000e-02,5.278000e-02,7.778000e-02,1.027800e-01
,1.277800e-01,1.527800e-01,1.777800e-01,2.027800e-01,2.277800e-01
,2.527800e-01,2.777800e-01,3.027800e-01,3.277800e-01,3.527800e-01]

y=[0.44236,0.4308,0.42299,0.41427,0.40548,0.39908,0.39039,0.3845,0.37882
,0.37411,0.36759,0.36434,0.35864,0.35508,0.35138]

xdata = np.array(x)
ydata = np.array(y)

parameters, covariance = scipy.optimize.curve_fit(Srt, xdata, ydata)

NOTE:
The np.abs inside the function does not make sense, but the complex result of lambertw apparently can be complex. In this case an error is raised as there is no safe casting rule, causing curvefit to abort.

Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

Mabe reasonable initial conditions for Smax, s0, Km, Vmax provided as keyword argument to curve_fit like e.g. p0=(1., 0., 10. ,10.) can fix the issue with the result of lambertw being complex
Thanks, you solved my problem. I will try different initial conditions with my "real" data - there might be some challenges in there, where I have to hack the fit ...
0

Your first error is produced by the t+t0 expression. It t is a list x, that's a list "concatenate" expression, which is fine for [1,2,3]+[4,5] but not [1,2,3]+5. That's why x and y have to arrays.

In the second error, what did the

print("s",type(s))
print("s",s)

show? Apparently s is not an array, or even a list.

Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.