CODEWITHSUNDEEP
pythoncurve-fittingnumerical-methodsnumerical-integration
Joshua
Joshua

Reputation: 1079

Python - How can I fit a curve to a function that contains a numerically calculated integral?

I have the following code:

import numpy as np
import scipy.integrate as spi
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import math as mh

def GUFunction(z, Omega_Lambda):
    integral = spi.quad(lambda zvar: AuxIntegrandum(zvar, Omega_Lambda), 0.0, z)[0]
    DL = (1+z) * c/H0 * integral *1000000
    return (5*(mh.log(DL,10)-1))

def AuxIntegrandum(z, Omega_Lambda):
    Omega_m = 1 - Omega_Lambda

    return 1 / mh.sqrt(Omega_m*(1+z)**3 + Omega_Lambda)

def DataFit(filename):
    print curve_fit(GUFunction, ComputeData(filename)[0], ComputeData(filename)[1])

DataFit("data.dat")

data.dat has z values in the first column and GUF(z) values in the second column.

Upon executing this code, the compiler tells me that comparing an array to a value (+inf or -inf) is ambiguous.
I think this refers to the integration boundaries, where it looks to see if I want to integrate to infinity. For some reason it apparently puts all z-values from the data file into the integration boundary.
Is there some trick I don't know about which allows you to fit a curve to a numerically integrated function?

Here's the exact error:

Traceback (most recent call last):
  File "plot.py", line 83, in <module>
    DataFit("data.dat")
  File "plot.py", line 67, in DataFit
    print curve_fit(GUFunction, ComputeData(filename)[0], ComputeData(filename)[1])
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 736, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 377, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 454, in func_wrapped
    return func(xdata, *params) - ydata
  File "plot.py", line 57, in GUFunction
    integral = spi.quad(lambda zvar: AuxIntegrandum(zvar, Omega_Lambda), 0.0, z)[0]
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 323, in quad
    points)
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 372, in _quad
    if (b != Inf and a != -Inf):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Upvotes: 0

Views: 879

Answers (1)

ev-br
ev-br

Reputation: 26040

Short answer: curve_fit tries to evaluate the target function on an array of xdata, but quad cannot accept a vector argument. You need to define your target function via e.g. a list comprehension over an input array.

Let's cook up a minimum reproducible example:

In [33]: xdata = np.linspace(0, 3, 11)

In [34]: ydata = xdata**3

In [35]: def integr(x):
    ...:     return quad(lambda t: t**2, 0, x)[0]
    ...: 

In [36]: def func(x, a):
    ...:     return integr(x) * a
    ...: 

In [37]: curve_fit(func, xdata, ydata)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-37-4660c65f85a2> in <module>()
----> 1 curve_fit(func, xdata, ydata)

 [... removed for clarity ...]

~/virtualenvs/py35/lib/python3.5/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    370 def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points):
    371     infbounds = 0
--> 372     if (b != Inf and a != -Inf):
    373         pass   # standard integration
    374     elif (b == Inf and a != -Inf):

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Which is precisely the error you're seeing. Ok, the error comes from quad, which tries to evaluate func(xdata, a), which boils down to integr(xdata) and that does not work. (How I found it out? I put import pdb; pdf.set_trace() inside the func function and poked around in the debugger).

Then, let's make the target function handle array arguments: 

In [38]: def func2(x, a):
    ...:     return np.asarray([integr(xx) for xx in x]) * a
    ...: 

In [39]: curve_fit(func2, xdata, ydata)
Out[39]: (array([ 3.]), array([[  3.44663413e-32]]))

Upvotes: 1

Related Questions