Optimize Function - use array as input

Question

I am playing with SciPy today and I wanted to test least square fitting. The function malo(time) works perfectly in returning me calculated concentrations if I put it in a loop which iterates over an array of timesteps (in the code "time").

Now I want to compare my calculated concentrations with my measured ones. I created a residuals function which calculates the difference between measured concentration (in the script an array called conc) and the modelled concentration with malo(time).

With optimize.leastsq I want to fit the parameter PD to fit both curves as good as possible. I don't see a mistake in my code, malo(time) performs well, but whenever I want to run the optimize.leastsq command Python says "only length-1 arrays can be converted to Python scalars". If I set the timedt array to a single value, the code runs without any error.

Do you see any chance to convince Python to use my array of timesteps in the loop?

import pylab as p
import math as m
import numpy as np
from scipy import optimize


Q = 0.02114
M = 7500.0
dt = 30.0
PD = 0.020242215
tom  = 26.0 #Minuten
tos = tom * 60.0 #Sekunden
timedt = np.array([30.,60.,90])
conc= np.array([ 2.7096,  2.258 ,  1.3548,  0.9032,  0.9032])


def malo(time):
     M1 = M/Q
     M2 = 1/(tos*m.sqrt(4*m.pi*PD*((time/tos)**3)))
     M3a = (1 - time/tos)**2
     M3b = 4*PD*(time/tos)
     M3 = m.exp(-1*(M3a/M3b))

     out = M1 * M2 * M3
     return out

def residuals(p,y,time):
    PD = p
    err = y - malo(timedt)
    return err

p0 = 0.05

p1 = optimize.leastsq(residuals,p0,args=(conc,timedt))

Answer 1

Notice that you're working here with arrays defined in NumPy module. Eg.

timedt = np.array([30.,60.,90])
conc= np.array([ 2.7096,  2.258 ,  1.3548,  0.9032,  0.9032])

Now, those arrays are not part of standard Python (which is a general purpose language). The problem is that you're mixing arrays with regular operations from the math module, which is part of the standard Python and only meant to work on scalars.

So, for example:

M2 = 1/(tos*m.sqrt(4*m.pi*PD*((time/tos)**3)))

will work if you use np.sqrt instead, which is designed to work on arrays:

M2 = 1/(tos*np.sqrt(4*m.pi*PD*((time/tos)**3)))

And so on.

NB: SciPy and other modules meant for numeric/scientific programming know about NumPy and are built on top of it, so those functions should all work on arrays. Just don't use math when working with them. NumPy comes with replicas of all those functions (sqrt, cos, exp, ...) to work with your arrays.