2

I am very new to Scipy and am already very impressed. I was hoping to use the ODEint and curve fitting to do some simple chemical kinetics modelling. Basically I am integrating an equation to get y values, then fitting this to some data that I am reading in from excel. For the simple process below it works perfect

This works:

data = pd.read_excel('C:\Spreadsheet', header = 0, parse_cols = 5)
data.head()
SM_CONC = 0.1176

def fitfuncSM(t, k1):
    def myode(SM, t):
        return -k1 * SM
    SM0 = SM_CONC
    SMt = odeint(myode, SM0, t)
    return SMt[:,0]
SM_fit, SMcov = curve_fit(fitfuncSM, data['Time'], data['A'], p0 = 1)

The problem

I want to add more complexity by curve fitting multiple plots at the same time (see below). This is a reaction of A converts to B converts to C. I am looking to curve fit and extract values k1 and k2. I have ydata for Ca and Cb. The function to be integrated works well (mytest - see below). However I am not sure how to import the ydata to be fitted. If I import an array of ydata's like this: SM_f, SMc = curve_fit(SM, data['Time'], ydata (array of the Ca and Cb), p0 = (1, 1). It gives me an error : operands could not be broadcast together with shapes (10,3) (3,10)"

def SM(t, k1, k2):
      def mytest(C, t):
        Ca = C[0]
        Cb = C[1]
        Cc = C[2]
        dAdt = -k1 * Ca
        dBdt = k1 * Ca - (k2 * Cb)
        dCdt = k2 * Cb
        dydt = [dAdt, dBdt, dCdt]
        return dydt
    init = [0.117, 0, 0]
    value = odeint(mytest, init, t)  
    return value

Update below - I have updated my code using the Transform you suggested (see below) I have the feeling I am missing something very simple. Any ideas? Is this a limitation of the curve fitting? Should I be using something else?

xdata = (data['Time'])
ydata = np.array([data['A'], data['B'], data['C']])
SM_CONC = 0.1176

def SM(t, k1, k2):
    def mytest(C, t):
        Ca = C[0]
        Cb = C[1]
        Cc = C[2]
        dAdt = -k1 * Ca
        dBdt = k1 * Ca - (k2 * Cb)
        dCdt = k2 * Cb
        dydt = np.array([dAdt, dBdt, dCdt])
        return dydt

    init = [SM_CONC, 0, 0]
    value = odeint(mytest, init, t) 
    return value[:,:]

SM_f, SMc = curve_fit(SM, xdata, ydata.T, p0 = (1, 1))

Error message

  ValueError                                Traceback (most recent call last)
  ValueError: object too deep for desired array

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-58-d09497ebe5ae> in <module>()
 30     return value[:,:]
 31 
---> 32 SM_f, SMc = curve_fit(SM, xdata, ydata.T, p0 = (1, 1))
 33 
 34 

C:\Users\Family\Anaconda3\lib\site-packages\scipy\optimize\minpack.py in curve_fit(f, xdata,          

553     # Remove full_output from kw, otherwise we're passing it in twice.
554     return_full = kw.pop('full_output', False)
--> 555     res = leastsq(func, p0, args=args, full_output=1, **kw)
556     (popt, pcov, infodict, errmsg, ier) = res
557 

C:\Users\Family\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)

377             maxfev = 200*(n + 1)
378         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 379                                  gtol, maxfev, epsfcn, factor, diag)
380     else:
381         if col_deriv:

error: Result from function call is not a proper array of floats.
DrPeter
  • 21
  • 3
  • Tips: use the code formatting option; with Python tracebacks, copy-paste the complete traceback (it's relatively clear here, but not always). –  Jan 10 '15 at 14:11
  • 1
    Probably `ydata` just needs to be transposed. Use `ydata.T` inside the `curve_fit(...)` call. –  Jan 10 '15 at 14:13
  • @Evert. I have updated with the code above. Any suggestions? – DrPeter Jan 11 '15 at 22:30
  • The error indicates `SM` should return an 1D array of floats; it looks like you're returning a 2D array. The problem is that your input `ydata` is multi-dimensional; only the inpudat `xdata` can be 2D. So you have to rephrase your problem. –  Jan 12 '15 at 07:40
  • Another way to look at the problem occurring is that your current `ydata` is `N x 3` dimensional, while the value you return from your fitfunction is `N x 2`; you can't subtract those two arrays, thus you can't calculate a chi-sq and the fit fails. –  Jan 12 '15 at 07:44
  • @evert I got around the problem by using Scipy Orthogonal distance regression instead and it worked as I wanted! I can post the code I used. Thanks very much for your help. – DrPeter Jan 21 '15 at 11:38
  • You're allowed (encouraged) to answer your own question if there's no appropriate answer, and you deem it may be useful to know for other people in the future. –  Jan 21 '15 at 12:18

0 Answers0